Geometric Deep Learning
Grids, Groups, Graphs,
Geodesics, and Gauges

Preface

For nearly two millenia since Euclid’s Elements, the word‘geometry’ has been synonymous with Euclidean geometry, as no other types of geometry existed. Euclid’s monopoly came to an end in the nineteenth century, with examples of non-Euclidean geometries constructed by Lobachevesky, Bolyai, Gauss, and Riemann. Towards the end of that century, these studies had diverged into disparate fields, with mathematicians and philosophers debating the validity of and relations between these geometries as well as the nature of the “one true geometry”.

A way out of this pickle was shown by a young mathematician Felix Klein, appointed in 1872 as professor in the small Bavarian University of Erlangen. In a research prospectus, which entered the annals of mathematics as the Erlangen Programme, Klein proposed approaching geometry as the study of invariants, i.e. properties unchanged under some class of transformations, called the symmetries of the geometry. This approach created clarity by showing that various geometries known at the time could be defined by an appropriate choice of symmetry transformations, formalized using the language of group theory. For instance, Euclidean geometry is concerned with lengths and angles, because these properties are preserved by the group of Euclidean transformations (rotations and translations), while affine geometry studies parallelism, which is preserved by the group of affine transformations. The relation between these geometries is immediately apparent when considering the respective groups, because the Euclidean group is a subgroup of the affine group, which in turn is a subgroup of the group of projective transformations.

The impact of the Erlangen Programme on geometry was very profound. Furthermore, it spilled to other fields, especially physics, where symmetry principles allowed to derive conservation laws from first principles of symmetry (an astonishing result known as Noether’s Theorem), and even enabled the classification of elementary particles as irreducible representations of the symmetry group. Category theory, now pervasive in pure mathematics, can be “regarded as a continuation of the Klein Erlangen Programme, in the sense that a geometrical space with its group of transformations is generalized to a category with its algebra of mappings”, in the words of its creators Samuel Eilenber and Saunders Mac Lane.

At the time of writing, the state of the field of deep learning is somewhat reminiscent of the field of geometry in the nineteenth century. There is a veritable zoo of neural network architectures for various kinds of data, but few unifying principles. As in times past, this makes it difficult to understand the relations between various methods, inevitably resulting in the reinvention and re-branding of the same concepts in different application domains. For a novice trying to learn the field, absorbing the sheer volume of redundant ideas is a true nightmare.

In this text, we make a modest attempt to apply the Erlangen Programme mindset to the domain of deep learning, with the ultimate goal of obtaining a systematisation of this field and ‘connecting the dots’. We call this geometrisation attempt ‘Geometric Deep Learning’, and true to the spirit of Felix Klein, propose to derive different inductive biases and network architectures implementing them from first principles of symmetry and invariance. In particular, we focus on a large class of neural networks designed for analysing unstructured sets, grids, graphs, and manifolds, and show that they can be understood in a unified manner as methods that respect the structure and symmetries of these domains.

We believe this text would appeal to a broad audience of deep learning researchers, practitioners, and enthusiasts. A novice may use it as an overview and introduction to Geometric Deep Learning. A seasoned deep learning expert may discover new ways of deriving familiar architectures from basic principles and perhaps some surprising connections. Practitioners may get new insights on how to solve problems in their respective fields.

With such a fast-paced field as modern machine learning, the risk of writing a text like this is that it becomes obsolete and irrelevant before it sees the light of day. Having focused on foundations, our hope is that the key concepts we discuss will transcend their specific realisations— or, as Claude Adrien Helvétius put it, “la connaissance de certains principes supplée facilement à la connoissance de certains faits.”

Notation

$\Omega,u$	Domain, point on domain
$x(u)\in\mathcal{X}(\Omega,\mathcal{C})$	Signal on the domain of the form $x:\Omega\rightarrow\mathcal{C}$
$f(x)\in\mathcal{F}(\mathcal{X}(\Omega))$	Functions on signals on the domain of the form $f:\mathcal{X}(\Omega)\rightarrow\mathcal{Y}$
${\mathfrak{G}},{\mathfrak{g}}$	Group, element of the group
${\mathfrak{g}}.u,\rho({\mathfrak{g}})$	Group action, group representation
$\mathbf{X}\in\mathcal{C}^{|\Omega|\times s}$	Matrix representing a signal on a discrete domain
$\mathbf{x}_{u}\in\mathcal{C}^{s}$	Vector representing a discrete domain signal $\mathbf{X}$ on element $u\in\Omega$
$x_{uj}\in\mathcal{C}$	Scalar representing the $j$th component of a discrete domain signal $\mathbf{X}$ on element $u\in\Omega$
$\mathbf{F}(\mathbf{X})$	Function on discrete domain signals that returns another discrete domain signal, as a matrix
$\tau:\Omega\rightarrow\Omega$	Automorphism of the domain
$\eta:\Omega\rightarrow\Omega^{\prime}$	Isomorphism between two different domains
$\sigma:\mathcal{C}\rightarrow\mathcal{C}^{\prime}$	Activation function (point-wise non-linearity)
$G=(\mathcal{V},\mathcal{E})$	Graph with nodes $\mathcal{V}$ and edges $\mathcal{E}$
$\mathcal{T}=(\mathcal{V},\mathcal{E},\mathcal{F})$	Mesh with nodes $\mathcal{V}$, edges $\mathcal{E}$, and faces $\mathcal{F}$
$x\star\theta$	Convolution with filter $\theta$
$S_{v}$	Shift operator
$\varphi_{i}$	Basis function
$T_{u}\Omega,T\Omega$	Tangent space at $u$, tangent bundle
$X\in T_{u}\Omega$	Tangent vector
$g_{u}(X,Y)=\langle X,Y\rangle_{u}$	Riemannian metric
$\ell(\gamma),\ell_{uv}$	Length of a curve $\gamma$, discrete metric on edge $(u,v)$

1 Introduction

The last decade has witnessed an experimental revolution in data science and machine learning, epitomised by deep learning methods. Indeed, many high-dimensional learning tasks previously thought to be beyond reach – such as computer vision, playing Go, or protein folding – are in fact feasible with appropriate computational scale. Remarkably, the essence of deep learning is built from two simple algorithmic principles: first, the notion of representation or feature learning, whereby adapted, often hierarchical, features capture the appropriate notion of regularity for each task, and second, learning by local gradient-descent, typically implemented as backpropagation.

While learning generic functions in high dimensions is a cursed estimation problem, most tasks of interest are not generic, and come with essential pre-defined regularities arising from the underlying low-dimensionality and structure of the physical world. This text is concerned with exposing these regularities through unified geometric principles that can be applied throughout a wide spectrum of applications.

Exploiting the known symmetries of a large system is a powerful and classical remedy against the curse of dimensionality, and forms the basis of most physical theories. Deep learning systems are no exception, and since the early days researchers have adapted neural networks to exploit the low-dimensional geometry arising from physical measurements, e.g. grids in images, sequences in time-series, or position and momentum in molecules, and their associated symmetries, such as translation or rotation. Throughout our exposition, we will describe these models, as well as many others, as natural instances of the same underlying principle of geometric regularity.

Such a ‘geometric unification’ endeavour in the spirit of the Erlangen Program serves a dual purpose: on one hand, it provides a common mathematical framework to study the most successful neural network architectures, such as CNNs, RNNs, GNNs, and Transformers. On the other, it gives a constructive procedure to incorporate prior physical knowledge into neural architectures and provide principled way to build future architectures yet to be invented.

Before proceeding, it is worth noting that our work concerns representation learning architectures and exploiting the symmetries of data therein. The many exciting pipelines where such representations may be used (such as self-supervised learning, generative modelling, or reinforcement learning) are not our central focus. Hence, we will not review in depth influential neural pipelines such as variational autoencoders (Kingma and Welling, 2013), generative adversarial networks (Goodfellow et al., 2014), normalising flows (Rezende and Mohamed, 2015), deep Q-networks (Mnih et al., 2015), proximal policy optimisation (Schulman et al., 2017), or deep mutual information maximisation (Hjelm et al., 2019). That being said, we believe that the principles we will focus on are of significant importance in all of these areas.

Further, while we have attempted to cast a reasonably wide net in order to illustrate the power of our geometric blueprint, our work does not attempt to accurately summarise the entire existing wealth of research on Geometric Deep Learning. Rather, we study several well-known architectures in-depth in order to demonstrate the principles and ground them in existing research, with the hope that we have left sufficient references for the reader to meaningfully apply these principles to any future geometric deep architecture they encounter or devise.

2 Learning in High Dimensions

Supervised machine learning, in its simplest formalisation, considers a set of $N$ observations ${\mathcal{D}}=\{(x_{i},y_{i})\}_{i=1}^{N}$ drawn i.i.d. from an underlying data distribution $P$ defined over ${\mathcal{X}}\times{\mathcal{Y}}$, where ${\mathcal{X}}$ and ${\mathcal{Y}}$ are respectively the data and the label domains. The defining feature in this setup is that ${\mathcal{X}}$ is a high-dimensional space: one typically assumes ${\mathcal{X}}=\mathbb{R}^{d}$ to be a Euclidean space of large dimension $d$.

Let us further assume that the labels $y$ are generated by an unknown function $f$, such that $y_{i}=f(x_{i})$, and the learning problem reduces to estimating the function $f$ using a parametrised function class ${\mathcal{F}}=\{f_{\bm{\theta}\in\Theta}\}$. Neural networks are a common realisation of such parametric function classes, in which case $\bm{\theta}\in\Theta$ corresponds to the network weights. In this idealised setup, there is no noise in the labels, and modern deep learning systems typically operate in the so-called interpolating regime, where the estimated $\tilde{f}\in{\mathcal{F}}$ satisfies $\tilde{f}(x_{i})=f(x_{i})$ for all $i=1,\ldots,N$. The performance of a learning algorithm is measured in terms of the expected performance on new samples drawn from ${P}$, using some loss $L(\cdot,\cdot)$

$${\mathcal{R}}(\tilde{f}):=\mathbb{E}_{P}\,\,L(\tilde{f}(x),f(x)),$$

with the squared-loss $L(y,y^{\prime})=\frac{1}{2}|y-y^{\prime}|^{2}$ being among the most commonly used ones.

A successful learning scheme thus needs to encode the appropriate notion of regularity or inductive bias for $f$, imposed through the construction of the function class $\mathcal{F}$ and the use of regularisation. We briefly introduce this concept in the following section.

2.1 Inductive Bias via Function Regularity

Modern machine learning operates with large, high-quality datasets, which, together with appropriate computational resources, motivate the design of rich function classes ${\mathcal{F}}$ with the capacity to interpolate such large data. This mindset plays well with neural networks, since even the simplest choices of architecture yields a dense class of functions.The capacity to approximate almost arbitrary functions is the subject of various Universal Approximation Theorems; several such results were proved and popularised in the 1990s by applied mathematicians and computer scientists (see e.g. Cybenko (1989); Hornik (1991); Barron (1993); Leshno et al. (1993); Maiorov (1999); Pinkus (1999)).

Figure 1: Multilayer Perceptrons (Rosenblatt, 1958), the simplest feed-forward neural networks, are universal approximators: with just one hidden layer, they can represent combinations of step functions, allowing to approximate any continuous function with arbitrary precision.

Universal Approximation, however, does not imply an absence of inductive bias. Given a hypothesis space ${\mathcal{F}}$ with universal approximation, we can define a complexity measure $c:{\mathcal{F}}\to\mathbb{R}_{+}$ and redefine our interpolation problem as

$$\tilde{f}\in\arg\min_{g\in{\mathcal{F}}}c(g)\quad\mathrm{s.t.}\quad g(x_{i})=f(x_{i})\quad\mathrm{for}\,\,\,i=1,\ldots,N,$$

i.e., we are looking for the most regular functions within our hypothesis class. For standard function spaces, this complexity measure can be defined as a norm,making ${\mathcal{F}}$ a Banach space and allowing to leverage a plethora of theoretical results in functional analysis. In low dimensions, splines are a workhorse for function approximation. They can be formulated as above, with a norm capturing the classical notion of smoothness, such as the squared-norm of second-derivatives $\int_{-\infty}^{+\infty}|f^{\prime\prime}(x)|^{2}\mathrm{d}x$ for cubic splines.

In the case of neural networks, the complexity measure $c$ can be expressed in terms of the network weights, i.e. $c(f_{\bm{\theta}})={c}(\bm{\theta})$. The $L_{2}$-norm of the network weights, known as weight decay, or the so-called path-norm (Neyshabur et al., 2015) are popular choices in deep learning literature. From a Bayesian perspective, such complexity measures can also be interpreted as the negative log of the prior for the function of interest. More generally, this complexity can be enforced explicitly by incorporating it into the empirical loss (resulting in the so-called Structural Risk Minimisation), or implicitly, as a result of a certain optimisation scheme. For example, it is well-known that gradient-descent on an under-determined least-squares objective will choose interpolating solutions with minimal $L_{2}$ norm. The extension of such implicit regularisation results to modern neural networks is the subject of current studies (see e.g. Blanc et al. (2020); Shamir and Vardi (2020); Razin and Cohen (2020); Gunasekar et al. (2017)). All in all, a natural question arises: how to define effective priors that capture the expected regularities and complexities of real-world prediction tasks?

2.2 The Curse of Dimensionality

While interpolation in low-dimensions (with $d=1,2$ or $3$) is a classic signal processing task with very precise mathematical control of estimation errors using increasingly sophisticated regularity classes (such as spline interpolants, wavelets, curvelets, or ridgelets), the situation for high-dimensional problems is entirely different.

In order to convey the essence of the idea, let us consider a classical notion of regularity that can be easily extended to high dimensions: 1-Lipschitz- functions $f:\mathcal{X}\to\mathbb{R}$, i.e. functions satisfying $|f(x)-f(x^{\prime})|\leq\|x-x^{\prime}\|$ for all $x,x^{\prime}\in\mathcal{X}$. This hypothesis only asks the target function to be locally smooth, i.e., if we perturb the input $x$ slightly (as measured by the norm $\|x-x^{\prime}\|$), the output $f(x)$ is not allowed to change much. If our only knowledge of the target function $f$ is that it is $1$-Lipschitz, how many observations do we expect to require to ensure that our estimate $\tilde{f}$ will be close to $f$? Figure 2 reveals that the general answer is necessarily exponential in the dimension $d$, signaling that the Lipschitz class grows ‘too quickly’ as the input dimension increases: in many applications with even modest dimension $d$, the number of samples would be bigger than the number of atoms in the universe. The situation is not better if one replaces the Lipschitz class by a global smoothness hypothesis, such as the Sobolev Class ${\mathcal{H}}^{s}(\Omega_{d})$. Indeed, classic results (Tsybakov, 2008) establish a minimax rate of approximation and learning for the Sobolev class of the order $\epsilon^{-d/s}$, showing that the extra smoothness assumptions on $f$ only improve the statistical picture when $s\propto d$, an unrealistic assumption in practice.

Figure 2: We consider a Lipschitz function $f(x)=\sum_{j=1}^{2^{d}}z_{j}\phi(x-x_{j})$ where $z_{j}=\pm 1$, $x_{j}\in\mathbb{R}^{d}$ is placed in each quadrant, and $\phi$ a locally supported Lipschitz ‘bump’. Unless we observe the function in most of the $2^{d}$ quadrants, we will incur in a constant error in predicting it. This simple geometric argument can be formalised through the notion of Maximum Discrepancy (von Luxburg and Bousquet, 2004), defined for the Lipschitz class as $\kappa(d)=\mathbb{E}_{x,x^{\prime}}\sup_{f\in\mathrm{Lip}(1)}\left|\frac{1}{N}\sum_{l}f(x_{l})-\frac{1}{N}\sum_{l}f(x^{\prime}_{l})\right|\simeq N^{-1/d}$, which measures the largest expected discrepancy between two independent $N$-sample expectations. Ensuring that $\kappa(d)\simeq\epsilon$ requires $N=\Theta(\epsilon^{-d})$; the corresponding sample $\{x_{l}\}_{l}$ defines an $\epsilon$-net of the domain. For a $d$-dimensional Euclidean domain of diameter $1$, its size grows exponentially as $\epsilon^{-d}$.

Fully-connected neural networks define function spaces that enable more flexible notions of regularity, obtained by considering complexity functions $c$ on their weights. In particular, by choosing a sparsity-promoting regularisation, they have the ability to break this curse of dimensionality (Bach, 2017). However, this comes at the expense of making strong assumptions on the nature of the target function $f$, such as that $f$ depends on a collection of low-dimensional projections of the input (see Figure 3). In most real-world applications (such as computer vision, speech analysis, physics, or chemistry), functions of interest tend to exhibits complex long-range correlations that cannot be expressed with low-dimensional projections (Figure 3), making this hypothesis unrealistic. It is thus necessary to define an alternative source of regularity, by exploiting the spatial structure of the physical domain and the geometric priors of $f$, as we describe in the next Section 3.

Figure 3: If the unknown function $f$ is presumed to be well approximated as $f(\mathbf{x})\approx g(\mathbf{A}\mathbf{x})$ for some unknown $\mathbf{A}\in\mathbb{R}^{k\times d}$ with $k\ll d$, then shallow neural networks can capture this inductive bias, see e.g. Bach (2017). In typical applications, such dependency on low-dimensional projections is unrealistic, as illustrated in this example: a low-pass filter projects the input images to a low-dimensional subspace; while it conveys most of the semantics, substantial information is lost.

3 Geometric Priors

Modern data analysis is synonymous with high-dimensional learning. While the simple arguments of Section 2.1 reveal the impossibility of learning from generic high-dimensional data as a result of the curse of dimensionality, there is hope for physically-structured data, where we can employ two fundamental principles: symmetry and scale separation. In the settings considered in this text, this additional structure will usually come from the structure of the domain underlying the input signals: we will assume that our machine learning system operates on signals (functions) on some domain $\Omega$. While in many cases linear combinations of points on $\Omega$ is not well-defined, we can linearly combine signals on it, i.e., the space of signals forms a vector space. Moreover, since we can define an inner product between signals, this space is a Hilbert space.

As a typical illustration, take $\Omega=\mathbb{Z}_{n}\times\mathbb{Z}_{n}$ to be a two-dimensional $n\times n$ grid, $x$ an RGB image (i.e. a signal $x:\Omega\rightarrow\mathbb{R}^{3}$), and $f$ a function (such as a single-layer Perceptron) operating on $3n^{2}$-dimensional inputs. As we will see in the following with greater detail, the domain $\Omega$ is usually endowed with certain geometric structure and symmetries. Scale separation results from our ability to preserve important characteristics of the signal when transferring it onto a coarser version of the domain (in our example, subsampling the image by coarsening the underlying grid).

We will show that both principles, to which we will generically refer as geometric priors, are prominent in most modern deep learning architectures. In the case of images considered above, geometric priors are built into Convolutional Neural Networks (CNNs) in the form of convolutional filters with shared weights (exploiting translational symmetry) and pooling (exploiting scale separation). Extending these ideas to other domains such as graphs and manifolds and showing how geometric priors emerge from fundamental principles is the main goal of Geometric Deep Learning and the leitmotif of our text.

3.1 Symmetries, Representations, and Invariance

Informally, a symmetry of an object or system is a transformation that leaves a certain property of said object or system unchanged or invariant. Such transformations may be either smooth, continuous, or discrete. Symmetries are ubiquitous in many machine learning tasks. For example, in computer vision the object category is unchanged by shifts, so shifts are symmetries in the problem of visual object classification. In computational chemistry, the task of predicting properties of molecules independently of their orientation in space requires rotational invariance. Discrete symmetries emerge naturally when describing particle systems where particles do not have canonical ordering and thus can be arbitrarily permuted, as well as in many dynamical systems, via the time-reversal symmetry (such as systems in detailed balance or the Newton’s second law of motion). As we will see in Section 4.1, permutation symmetries are also central to the analysis of graph-structured data.

Symmetry groups

The set of symmetries of an object satisfies a number of properties. First, symmetries may be combined to obtain new symmetries: if ${\mathfrak{g}}$ and ${\mathfrak{h}}$ are two symmetries, then their compositions ${\mathfrak{g}}\circ{\mathfrak{h}}$ and ${\mathfrak{h}}\circ{\mathfrak{g}}$ are also symmetries. The reason is that if both transformations leave the object invariant, then so does the composition of transformations, and hence the composition is also a symmetry. Furthermore, symmetries are always invertible, and the inverse is also a symmetry. This shows that the collection of all symmetries form an algebraic object known as a group. Since these objects will be a centerpiece of the mathematical model of Geometric Deep Learning, they deserve a formal definition and detailed discussion:

Note that commutativity is not part of this definition, i.e. we may have ${\mathfrak{g}}{\mathfrak{h}}\neq{\mathfrak{h}}{\mathfrak{g}}$. Groups for which ${\mathfrak{g}}{\mathfrak{h}}={\mathfrak{h}}{\mathfrak{g}}$ for all ${\mathfrak{g}},{\mathfrak{h}}\in{\mathfrak{G}}$ are called commutative or Abelian.

Though some groups can be very large and even infinite, they often arise from compositions of just a few elements, called group generators. Formally, $\mathfrak{G}$ is said to be generated by a subset $S\subseteq\mathfrak{G}$ (called the group generator) if every element ${\mathfrak{g}}\in{\mathfrak{G}}$ can be written as a finite composition of the elements of $S$ and their inverses. For instance, the symmetry group of an equilateral triangle (dihedral group $\mathrm{D}_{3}$) is generated by a $60^{\circ}$ rotation and a reflection (Figure 4). The 1D translation group, which we will discuss in detail in the following, is generated by infinitesimal displacements; this is an example of a Lie group of differentiable symmetries.

Note that here we have defined a group as an abstract object, without saying what the group elements are (e.g. transformations of some domain), only how they compose. Hence, very different kinds of objects may have the same symmetry group. For instance, the aforementioned group of rotational and reflection symmetries of a triangle is the same as the group of permutations of a sequence of three elements (we can permute the corners in the triangle in any way using a rotation and reflection – see Figure 4).

  

Figure 4: Left: an equilateral triangle with corners labelled by $1,2,3$, and all possible rotations and reflections of the triangle. The group $\mathrm{D}_{3}$ of rotation/reflection symmetries of the triangle is generated by only two elements (rotation by $60^{\circ}$ R and reflection F) and is the same as the group $\Sigma_{3}$ of permutations of three elements. Right: the multiplication table of the group $\mathrm{D}_{3}$. The element in the row ${\mathfrak{g}}$ and column ${\mathfrak{h}}$ corresponds to the element ${\mathfrak{g}}{\mathfrak{h}}$.

Group Actions and Group Representations

Rather than considering groups as abstract entities, we are mostly interested in how groups act on data. Since we assumed that there is some domain $\Omega$ underlying our data, we will study how the group acts on $\Omega$ (e.g. translation of points of the plane), and from there obtain actions of the same group on the space of signals $\mathcal{X}(\Omega)$ (e.g. translations of planar images and feature maps).

A group action of ${\mathfrak{G}}$ on a set $\Omega$ is defined as a mapping $({\mathfrak{g}},u)\mapsto{\mathfrak{g}}.u$ associating a group element ${\mathfrak{g}}\in{\mathfrak{G}}$ and a point $u\in\Omega$ with some other point on $\Omega$ in a way that is compatible with the group operations, i.e., ${\mathfrak{g}}.({\mathfrak{h}}.u)=({\mathfrak{g}}{\mathfrak{h}}).u$ for all ${\mathfrak{g}},{\mathfrak{h}}\in{\mathfrak{G}}$ and $u\in\Omega$. We shall see numerous instances of group actions in the following sections. For example, in the plane the Euclidean group $\mathrm{E}(2)$ is the group of transformations of $\mathbb{R}^{2}$ that preserves Euclidean distances, and consists of translations, rotations, and reflections. The same group, however, can also act on the space of images on the plane (by translating, rotating and flipping the grid of pixels), as well as on the representation spaces learned by a neural network. More precisely, if we have a group ${\mathfrak{G}}$ acting on $\Omega$, we automatically obtain an action of ${\mathfrak{G}}$ on the space $\mathcal{X}(\Omega)$:

$$({\mathfrak{g}}.x)(u)=x({\mathfrak{g}}^{-1}u).$$ (3)

Due to the inverse on ${\mathfrak{g}}$, this is indeed a valid group action, in that we have $({\mathfrak{g}}.({\mathfrak{h}}.x))(u)=(({\mathfrak{g}}{\mathfrak{h}}).x)(u)$.

The most important kind of group actions, which we will encounter repeatedly throughout this text, are linear group actions, also known as group representations. The action on signals in equation (3) is indeed linear, in the sense that

$${\mathfrak{g}}.(\alpha x+\beta x^{\prime})=\alpha({\mathfrak{g}}.x)+\beta({\mathfrak{g}}.x^{\prime})$$

for any scalars $\alpha,\beta$ and signals $x,x^{\prime}\in\mathcal{X}(\Omega)$. We can describe linear actions either as maps $({\mathfrak{g}},x)\mapsto{\mathfrak{g}}.x$ that are linear in $x$, or equivalently, by currying, as a map $\rho:{\mathfrak{G}}\rightarrow\mathbb{R}^{n\times n}$that assigns to each group element ${\mathfrak{g}}$ an (invertible) matrix $\rho({\mathfrak{g}})$. The dimension $n$ of the matrix is in general arbitrary and not necessarily related to the dimensionality of the group or the dimensionality of $\Omega$, but in applications to deep learning $n$ will usually be the dimensionality of the feature space on which the group acts. For instance, we may have the group of 2D translations acting on a space of images with $n$ pixels.

As with a general group action, the assignment of matrices to group elements should be compatible with the group action. More specifically, the matrix representing a composite group element ${\mathfrak{g}}{\mathfrak{h}}$ should equal the matrix product of the representation of ${\mathfrak{g}}$ and ${\mathfrak{h}}$:

Written in the language of group representations, the action of ${\mathfrak{G}}$ on signals $x\in\mathcal{X}(\Omega)$ is defined as $\rho({\mathfrak{g}})x(u)=x({\mathfrak{g}}^{-1}u)$. We again verify that

$$(\rho({\mathfrak{g}})(\rho({\mathfrak{h}})x))(u)=(\rho({\mathfrak{g}}{\mathfrak{h}})x)(u).$$

Figure 5: Three spaces of interest in Geometric Deep Learning: the (physical) domain $\Omega$, the space of signals $\mathcal{X}(\Omega)$, and the hypothesis class $\mathcal{F}(\mathcal{X}(\Omega))$. Symmetries of the domain $\Omega$ (captured by the group ${\mathfrak{G}}$) act on signals $x\in\mathcal{X}(\Omega)$ through group representations $\rho({\mathfrak{g}})$, imposing structure on the functions $f\in\mathcal{F}(\mathcal{X}(\Omega))$ acting on such signals.

Invariant and Equivariant functions

The symmetry of the domain $\Omega$ underlying the signals $\mathcal{X}(\Omega)$ imposes structure on the function $f$ defined on such signals. It turns out to be a powerful inductive bias, improving learningefficiency by reducing the space of possible interpolants, $\mathcal{F}(\mathcal{X}(\Omega))$, to those which satisfy the symmetry priors. Two important cases we will be exploring in this text are invariant and equivariant functions.

A classical example of invariance is shift-invariance,arising in computer vision and pattern recognition applications such as image classification. The function $f$ in this case (typically implemented as a Convolutional Neural Network) inputs an image and outputs the probability of the image to contain an object from a certain class (e.g. cat or dog). It is often reasonably assumed that the classification result should not be affected by the position of the object in the image, i.e., the function $f$ must be shift-invariant. Multi-layer Perceptrons, which can approximate any smooth function, do not have this property – one of the reasons why early attempts to apply these architectures to problems of pattern recognition in the 1970s failed. The development of neural network architectures with local weight sharing, as epitomised by Convolutional Neural Networks, was, among other reasons, motivated by the need for shift-invariant object classification.

If we however take a closer look at the convolutional layers of CNNs, we will find that they are not shift-invariant but shift-equivariant: in other words, a shift of the input to a convolutional layer produces a shift in the output feature maps by the same amount.

Resorting again to computer vision, a prototypical application requiring shift-equivariance is image segmentation, where the output of $f$ is a pixel-wise image mask. Obviously, the segmentation mask must follow shifts in the input image. In this example, the domains of the input and output are the same, but since the input has three color channels while the output has one channel per class, the representations $(\rho,\mathcal{X}(\Omega,\mathcal{C}))$ and $(\rho^{\prime},\mathcal{X}(\Omega,\mathcal{C}^{\prime}))$ are somewhat different.

However, even the previous use case of image classification is usually implemented as a sequence of convolutional (shift-equivariant) layers, followed by global pooling (which is shift-invariant). As we will see in Section 3.5, this is a general blueprint of a majority of deep learning architectures, including CNNs and Graph Neural Networks (GNNs).

3.2 Isomorphisms and Automorphisms

Subgroups and Levels of structure

As mentioned before, a symmetryis a transformation that preserves some property or structure, and the set of all such transformations for a given structure forms a symmetry group. It happens often that there is not one but multiple structures of interest, and so we can consider several levels of structure on our domain $\Omega$. Hence, what counts as a symmetry depends on the structure under consideration, but in all cases a symmetry is an invertible map that respects this structure.

On the most basic level, the domain $\Omega$ is a set, which has a minimal amount of structure: all we can say is that the set has some cardinality. Self-maps that preserve this structure are bijections (invertible maps), which we may consider as set-level symmetries. One can easily verify that this is a group by checking the axioms: a compositions of two bijections is also a bijection (closure), the associativity stems from the associativity of the function composition, the map $\tau(u)=u$ is the identity element, and for every $\tau$ the inverse exists by definition, satisfying $(\tau\circ\tau^{-1})(u)=(\tau^{-1}\circ\tau)(u)=u$.

Depending on the application, there may be further levels of structure. For instance, if $\Omega$ is a topological space, we can consider maps that preserve continuity: such maps are called homeomorphisms and in addition to simple bijections between sets, are also continuous and have continuous inverse. Intuitively, continuous functions are well-behaved and map points in a neighbourhood (open set) around a point $u$ to a neighbourhood around $\tau(u)$.

One can further demand that the map and its inverse are (continuously) differentiable,i.e., the map and its inverse have a derivative at every point (and the derivative is also continuous). This requires further differentiable structure that comes with differentiable manifolds, where such maps are called diffeomorphisms and denoted by $\mathrm{Diff}(\Omega)$. Additional examples of structures we will encounter include distances or metrics (maps preserving them are called isometries) or orientation (to the best of our knowledge, orientation-preserving maps do not have a common Greek name).

The right level of structure to consider depends on the problem. For example, when segmenting histopathology slide images, we may wish to consider flipped versions of an image as equivalent (as the sample can be flipped when put under the microscope), but if we are trying to classify road signs, we would only want to consider orientation-preserving transformations as symmetries (since reflections could change the meaning of the sign).

As we add levels of structure to be preserved, the symmetry group will get smaller. Indeed, adding structure is equivalent to selecting a subgroup, which is a subset of the larger group that satisfies the axioms of a group by itself:

For instance, the group of Euclidean isometries $\operatorname{E}(2)$ is a subgroup of the group of planar diffeomorphisms $\operatorname{Diff}(2)$, and in turn the group of orientation-preserving isometries $\operatorname{SE}(2)$ is a subgroup of $\operatorname{E}(2)$. This hierarchy of structure follows the Erlangen Programme philosophy outlined in the Preface: in Klein’s construction, the Projective, Affine, and Euclidean geometries have increasingly more invariants and correspond to progressively smaller groups.

Isomorphisms and Automorphisms

We have described symmetries as structure preserving and invertible maps from an object to itself. Such maps are also known as automorphisms, and describe a way in which an object is equivalent it itself. However, an equally important class of maps are the so-called isomorphisms, which exhibit an equivalence between two non-identical objects. These concepts are often conflated, but distinguishing them is necessary to create clarity for our following discussion.

To understand the difference, consider a set $\Omega=\{0,1,2\}$. An automorphism of the set $\Omega$ is a bijection $\tau:\Omega\rightarrow\Omega$ such as a cyclic shift $\tau(u)=u+1\mod 3$. Such a map preserves the cardinality property, and maps $\Omega$ onto itself. If we have another set $\Omega^{\prime}=\{a,b,c\}$ with the same number of elements, then a bijection $\eta:\Omega\rightarrow\Omega^{\prime}$ such as $\eta(0)=a$, $\eta(1)=b$, $\eta(2)=c$ is a set isomorphism.

As we will see in Section 4.1 for graphs, the notion of structure includes not just the number of nodes, but also the connectivity. An isomorphism $\eta:\mathcal{V}\rightarrow\mathcal{V}^{\prime}$ between two graphs $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and $\mathcal{G}^{\prime}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime})$ is thus a bijection between the nodes that maps pairs of connected nodes to pairs of connected nodes, and likewise for pairs of non-connected nodes.Two isomorphic graphs are thus structurally identical, and differ only in the way their nodes are ordered. On the other hand, a graph automorphism or symmetry is a map $\tau:\mathcal{V}\rightarrow\mathcal{V}$ maps the nodes of the graph back to itself, while preserving the connectivity. A graph with a non-trivial automorphism (i.e., $\tau\neq\mathrm{id}$) presents symmetries.

3.3 Deformation Stability

The symmetry formalism introduced in Sections 3.1–3.2 captures an idealised world where we know exactly which transformations are to be considered as symmetries, and we want to respect these symmetries exactly. For instance in computer vision, we might assume that planar translations are exact symmetries. However, the real world is noisy and this model falls short in two ways.

Firstly, while these simple groups provide a way to understand global symmetries of the domain $\Omega$ (and by extension, of signals on it, ${\mathcal{X}}(\Omega)$), they do not capture local symmetries well. For instance, consider a video scene with several objects, each moving along its own different direction. At subsequent frames, the resulting scene will contain approximately the same semantic information, yet no global translation explains the transformation from one frame to another. In other cases, such as a deformable 3D object viewed by a camera, it is simply very hard to describe the group of transformations that preserve the object identity. These examples illustrate that in reality we are more interested in a far larger set of transformations where global, exact invariance is replaced by a local, inexact one. In our discussion, we will distinguish between two scenarios: the setting where the domain $\Omega$ is fixed, and signals $x\in{\mathcal{X}}(\Omega)$ are undergoing deformations, and the setting where the domain $\Omega$ itself may be deformed.

Stability to signal deformations

In many applications, we know a priori that a small deformation of the signal $x$ should not change the output of $f(x)$, so it is tempting to consider such deformations as symmetries. For instance, we could view small diffeomorphisms $\tau\in\operatorname{Diff}(\Omega)$, or even small bijections, as symmetries. However, small deformations can be composed to form large deformations, so “small deformations” do not form a group,and we cannot ask for invariance or equivariance to small deformations only. Since large deformations can can actually materially change the semantic content of the input, it is not a good idea to use the full group $\operatorname{Diff}(\Omega)$ as symmetry group either.

A better approach is to quantify how “far” a given $\tau\in\operatorname{Diff}(\Omega)$ is from a given symmetry subgroup ${\mathfrak{G}}\subset\operatorname{Diff}(\Omega)$ (e.g. translations) with a complexity measure $c(\tau)$, so that $c(\tau)=0$ whenever $\tau\in{\mathfrak{G}}$. We can now replace our previous definition of exact invariance and equivarance under group actions with a ‘softer’ notion of deformation stability (or approximate invariance):

$$\|f(\rho(\tau)x)-f(x)\|\leq Cc(\tau)\|x\|,~{},~{}\forall x\in{\mathcal{X}}(\Omega)$$ (4)

where $\rho(\tau)x(u)=x(\tau^{-1}u)$ as before, and where $C$ is some constant independent of the signal $x$. A function $f\in\mathcal{F}(\mathcal{X}(\Omega))$ satisfying the above equation is said to be geometrically stable. We will see examples of such functions in the next Section 3.4.

Since $c(\tau)=0$ for $\tau\in{\mathfrak{G}}$, this definition generalises the ${\mathfrak{G}}$-invariance property defined above. Its utility in applications depends on introducing an appropriate deformation cost. In the case of images defined over a continuous Euclidean plane, a popular choice is $c^{2}(\tau):=\int_{\Omega}\|\nabla\tau(u)\|^{2}\mathrm{d}u$, which measures the ‘elasticity’ of $\tau$, i.e., how different it is from the displacement by a constant vector field. This deformation cost is in fact a norm often called the Dirichlet energy, and can be used to quantify how far $\tau$ is from the translation group.

Figure 6: The set of all bijective mappings from $\Omega$ into itself forms the set automorphism group $\mathrm{Aut}(\Omega)$, of which a symmetry group ${\mathfrak{G}}$ (shown as a circle) is a subgroup. Geometric Stability extends the notion of ${\mathfrak{G}}$-invariance and equivariance to ‘transformations around ${\mathfrak{G}}$’ (shown as gray ring), quantified in the sense of some metric between transformations. In this example, a smooth distortion of the image is close to a shift.

Stability to domain deformations

In many applications, the object being deformed is not the signal, but the geometric domain $\Omega$ itself. Canonical instances of this are applications dealing with graphs and manifolds: a graph can model a social network at different instance of time containing slightly different social relations (follow graph), or a manifold can model a 3D object undergoing non-rigid deformations. This deformation can be quantified as follows. If $\mathcal{D}$ denotes the space of all possible variable domains (such as the space of all graphs, or the space of Riemannian manifolds), one can define for $\Omega,\tilde{\Omega}\in\mathcal{D}$ an appropriate metric (‘distance’) $d(\Omega,\tilde{\Omega})$ satisfying $d(\Omega,\tilde{\Omega})=0$ if $\Omega$ and $\tilde{\Omega}$ are equivalent in some sense: for example, the graph edit distance vanishes when the graphs are isomorphic, and the Gromov-Hausdorff distance between Riemannian manifolds equipped with geodesic distances vanishes when two manifolds are isometric.

A common construction of such distances between domains relies on some family of invertible mapping $\eta:\Omega\to\tilde{\Omega}$ that try to ‘align’ the domains in a way that the corresponding structures are best preserved. For example, in the case of graphs or Riemannian manifolds (regarded as metric spaces with the geodesic distance), this alignment can compare pair-wise adjacency or distance structures ($d$ and $\tilde{d}$, respectively),

$$d_{{\mathcal{D}}}(\Omega,\tilde{\Omega})=\inf_{\eta\in{\mathfrak{G}}}\|d-\tilde{d}\circ(\eta\times\eta)\|$$

where ${\mathfrak{G}}$ is the group of isomorphisms such as bijections or isometries, and the norm is defined over the product space $\Omega\times\Omega$. In other words, a distance between elements of $\Omega,\tilde{\Omega}$ is ‘lifted’ to a distance between the domains themselves, by accounting for all the possible alignments that preserve the internal structure. Given a signal $x\in{\mathcal{X}}(\Omega)$ and a deformed domain $\tilde{\Omega}$, one can then consider the deformed signal $\tilde{x}=x\circ\eta^{-1}\in{\mathcal{X}}(\tilde{\Omega})$.

By slightly abusing the notation, we define ${\mathcal{X}}(\mathcal{D})=\{({\mathcal{X}}(\Omega),\Omega)\,:\,\Omega\in\mathcal{D}\}$ as the ensemble of possible input signals defined over a varying domain. A function $f:{\mathcal{X}}(\mathcal{D})\to{\mathcal{Y}}$ is stable to domain deformations if

$$\|f(x,\Omega)-f(\tilde{x},\tilde{\Omega})\|\leq C\|x\|d_{{\mathcal{D}}}(\Omega,\tilde{\Omega})~{}$$ (5)

for all $\Omega,\tilde{\Omega}\in\mathcal{D}$, and $x\in\mathcal{X}(\Omega)$. We will discuss this notion of stability in the context of manifolds in Sections 4.4–4.6, where isometric deformations play a crucial role. Furthermore, it can be shown that the stability to domain deformations is a natural generalisation of the stability to signal deformations, by viewing the latter in terms of deformations of the volume form Gama et al. (2019).

3.4 Scale Separation

While deformation stability substantially strengthens the global symmetry priors, it is not sufficient in itself to overcome the curse of dimensionality, in the sense that, informally speaking, there are still “too many" functions that respect (4) as the size of the domain grows. A key insight to overcome this curse is to exploit the multiscale structure of physical tasks. Before describing multiscale representations, we need to introduce the main elements of Fourier transforms, which rely on frequency rather than scale.

Fourier Transform and Global invariants

Arguably the most famous signal decomposition is the Fourier transform, the cornerstone of harmonic analysis. The classical one-dimensional Fourier transform

$$\hat{x}(\xi)=\int_{-\infty}^{+\infty}x(u)e^{-\mathrm{i}\xi u}\mathrm{d}u$$

expresses the function $x(u)\in L^{2}(\Omega)$ on the domain $\Omega=\mathbb{R}$ as a linear combination of orthogonal oscillating basis functions $\varphi_{\xi}(u)=e^{\mathrm{i}\xi u}$, indexed by their rate of oscillation (or frequency) $\xi$. Such an organisation into frequencies reveals important information about the signal, e.g. its smoothness and localisation. The Fourier basis itself has a deep geometric foundation and can be interpreted as the natural vibrations of the domain, related to its geometric structure (see e.g. Berger (2012)).

The Fourier transformplays a crucial role in signal processing as it offers a dual formulation of convolution,

$$(x\star\theta)(u)=\int_{-\infty}^{+\infty}x(v)\theta(u-v)\mathrm{d}v$$

a standard model of linear signal filtering (here and in the following, $x$ denotes the signal and $\theta$ the filter). As we will show in the following, the convolution operator is diagonalised in the Fourier basis, making it possible to express convolution as the product of the respective Fourier transforms,

$$\widehat{(x\star\theta)}(\xi)=\hat{x}(\xi)\cdot\hat{\theta}(\xi),$$

a fact known in signal processing as the Convolution Theorem.

As it turns out, many fundamental differential operators such as the Laplacian are described as convolutions on Euclidean domains. Since such differential operators can be defined intrinsically over very general geometries, this provides a formal procedure to extend Fourier transforms beyond Euclidean domains, including graphs, groups and manifolds. We will discuss this in detail in Section 4.4.

An essential aspect of Fourier transforms is that they reveal global properties of the signal and the domain, such as smoothness or conductance. Such global behavior is convenient in presence of global symmetries of the domain such as translation, but not to study more general diffeomorphisms. This requires a representation that trades off spatial and frequential localisation, as we see next.

Multiscale representations

The notion of local invariance can be articulated by switching from a Fourier frequency-based representation to a scale-based representation, the cornerstone of multi-scale decomposition methods such as wavelets.The essential insight of multi-scale methods is to decompose functions defined over the domain $\Omega$ into elementary functions that are localised both in space and frequency.In the case of wavelets, this is achieved by correlating a translated and dilated filter (mother wavelet) $\psi$, producing a combined spatio-frequency representation called a continuous wavelet transform

$$(W_{\psi}x)(u,\xi)=\xi^{-1/2}\int_{-\infty}^{+\infty}\psi\left(\frac{v-u}{\xi}\right)x(v)\mathrm{d}v.$$

The translated and dilated filters are called wavelet atoms; their spatial position and dilation correspond to the coordinates $u$ and $\xi$ of the wavelet transform. These coordinates are usually sampled dyadically ($\xi=2^{-j}$ and $u=2^{-j}k$), with $j$ referred to as scale. Multi-scale signal representations bring important benefits in terms of capturing regularity properties beyond global smoothness, such as piece-wise smoothness, which made them a popular tool in signal and image processing and numerical analysis in the 90s.

Deformation stability of Multiscale representations:

The benefit of multiscale localised wavelet decompositions over Fourier decompositions is revealed when considering the effect of small deformations ‘nearby’ the underlying symmetry group. Let us illustrate this important concept in the Euclidean domain and the translation group. Since the Fourier representation diagonalises the shift operator (which can be thought of as convolution, as we will see in more detail in Section 4.2), it is an efficient representation for translation transformations. However, Fourier decompositions are unstable under high-frequency deformations. In contrast, wavelet decompositions offer a stable representation in such cases.

Indeed, let us consider $\tau\in\mathrm{Aut}(\Omega)$ and its associated linear representation $\rho(\tau)$. When $\tau(u)=u-v$ is a shift, as we will verify in Section 4.2, the operator $\rho(\tau)=S_{v}$ is a shift operator that commutes with convolution. Since convolution operators are diagonalised by the Fourier transform, the action of shift in the frequency domain amounts to shifting the complex phase of the Fourier transform,

$$(\widehat{S_{v}x})(\xi)=e^{-\mathrm{i}\xi v}\hat{x}(\xi).$$

Thus, the Fourier modulus $f(x)=|\hat{x}|$ removing the complex phase is a simple shift-invariant function, $f(S_{v}x)=f(x)$. However, if we have only approximate translation, $\tau(u)=u-\tilde{\tau}(u)$ with $\|\nabla\tau\|_{\infty}=\sup_{u\in\Omega}\|\nabla\tilde{\tau}(u)\|\leq\epsilon$, the situation is entirely different: it is possible to show that

$$\frac{\|f(\rho(\tau)x)-f(x)\|}{\|x\|}=\mathcal{O}(1)$$

irrespective of how small $\epsilon$ is (i.e., how close is $\tau$ to being a shift). Consequently, such Fourier representation is unstable under deformations, however small. This unstability is manifested in general domains and non-rigid transformations; we will see another instance of this unstability in the analysis of 3d shapes using the natural extension of Fourier transforms described in Section 4.4.

Wavelets offer a remedy to this problem that also reveals the power of multi-scale representations. In the above example, we can show (Mallat, 2012) that the wavelet decomposition $W_{\psi}x$ is approximately equivariant to deformations,

$$\frac{\|\rho(\tau)(W_{\psi}x)-W_{\psi}(\rho(\tau)x)\|}{\|x\|}=\mathcal{O}(\epsilon).$$

In other words, decomposing the signal information into scales using localised filters rather than frequencies turns a global unstable representation into a family of locally stable features. Importantly, such measurements at different scales are not yet invariant, and need to be progressively processed towards the low frequencies, hinting at the deep compositional nature of modern neural networks, and captured in our Blueprint for Geometric Deep Learning, presented next.

Scale Separation Prior:

We can build from this insight by considering a multiscale coarsening of the data domain $\Omega$ into a hierarchy $\Omega_{1},\ldots,\Omega_{J}$. As it turns out, such coarsening can be defined on very general domains, including grids, graphs, and manifolds. Informally, a coarsening assimilates nearby points $u,u^{\prime}\in\Omega$ together, and thus only requires an appropriate notion of metric in the domain. If ${\mathcal{X}}_{j}(\Omega_{j},\mathcal{C}_{j}):=\{x_{j}:\Omega_{j}\to\mathcal{C}_{j}\}$ denotes signals defined over the coarsened domain $\Omega_{j}$, we informally say that a function $f:{\mathcal{X}}(\Omega)\to{\mathcal{Y}}$ is locally stable at scale $j$ if it admits a factorisation of the form $f\approx f_{j}\circ P_{j}$, where $P_{j}:{\mathcal{X}}(\Omega)\to{\mathcal{X}}_{j}(\Omega_{j})$ is a non-linear coarse graining and $f_{j}:{\mathcal{X}}_{j}(\Omega_{j})\to{\mathcal{Y}}$. In other words, while the target function $f$ might depend on complex long-range interactions between features over the whole domain, in locally-stable functions it is possible to separate the interactions across scales, by first focusing on localised interactions that are then propagated towards the coarse scales.

Figure 7: Illustration of Scale Separation for image classification tasks. The classifier $f^{\prime}$ defined on signals on the coarse grid $\mathcal{X}(\Omega^{\prime})$ should satisfy $f\approx f^{\prime}\circ P$, where $P:\mathcal{X}(\Omega)\rightarrow\mathcal{X}(\Omega^{\prime})$.

Such principlesare of fundamental importance in many areas of physics and mathematics, as manifested for instance in statistical physics in the so-called renormalisation group, or leveraged in important numerical algorithms such as the Fast Multipole Method. In machine learning, multiscale representations and local invariance are the fundamental mathematical principles underpinning the efficiency of Convolutional Neural Networks and Graph Neural Networks and are typically implemented in the form of local pooling. In future work, we will further develop tools from computational harmonic analysis that unify these principles across our geometric domains and will shed light onto the statistical learning benefits of scale separation.

3.5 The Blueprint of Geometric Deep Learning

The geometric principles of Symmetry, Geometric Stability, and Scale Separation discussed in Sections 3.1–3.4 can be combined to provide a universal blueprint for learning stable representations of high-dimensional data. These representations will be produced by functions $f$ operating on signals $\mathcal{X}(\Omega,\mathcal{C})$ defined on the domain $\Omega$, which is endowed with a symmetry group ${\mathfrak{G}}$.

The geometric priors we have described so far do not prescribe a specific architecture for building such representation, but rather a series of necessary conditions. However, they hint at an axiomatic construction that provably satisfies these geometric priors, while ensuring a highly expressive representation that can approximate any target function satisfying such priors.

A simple initial observation is that, in order to obtain a highly expressive representation, we are required to introduce a non-linear element, since if $f$ is linear and ${\mathfrak{G}}$-invariant, then for all $x\in{\mathcal{X}}(\Omega)$,

$$f(x)=\frac{1}{\mu({\mathfrak{G}})}\int_{{\mathfrak{G}}}f({\mathfrak{g}}.x)\mathrm{d}\mu({\mathfrak{g}})=f\left(\frac{1}{\mu({\mathfrak{G}})}\int_{{\mathfrak{G}}}({\mathfrak{g}}.x)\mathrm{d}\mu({\mathfrak{g}})\right),$$

which indicates that $F$ only depends on $x$ through the ${\mathfrak{G}}$-average $A{x}=\frac{1}{\mu({\mathfrak{G}})}\int_{{\mathfrak{G}}}({\mathfrak{g}}.x)\mathrm{d}\mu({\mathfrak{g}})$. In the case of images and translation, this would entail using only the average RGB color of the input!

While this reasoning shows that the family of linear invariants is not a very rich object, the family of linear equivariants provides a much more powerful tool, since it enables the construction of rich and stable features by composition with appropriate non-linear maps, as we will now explain. Indeed, if $B:{\mathcal{X}}(\Omega,{\mathcal{C}})\to{\mathcal{X}}(\Omega,{\mathcal{C}}^{\prime})$ is ${\mathfrak{G}}$-equivariant satisfying $B({\mathfrak{g}}.x)={\mathfrak{g}}.B(x)$ for all $x\in{\mathcal{X}}$ and ${\mathfrak{g}}\in{\mathfrak{G}}$, and $\sigma:{\mathcal{C}}^{\prime}\to{\mathcal{C}}^{\prime\prime}$ is an arbitrary (non-linear) map, then we easily verify that the composition $U:=(\bm{\sigma}\circ B):{\mathcal{X}}(\Omega,{\mathcal{C}})\to{\mathcal{X}}(\Omega,{\mathcal{C}}^{\prime\prime})$ is also ${\mathfrak{G}}$-equivariant, where $\bm{\sigma}:{\mathcal{X}}(\Omega,{\mathcal{C}}^{\prime})\to{\mathcal{X}}(\Omega,{\mathcal{C}}^{\prime\prime})$ is the element-wise instantiation of $\sigma$ given as $(\bm{\sigma}(x))(u):=\sigma(x(u))$.

This simple property allows us to define a very general family of ${\mathfrak{G}}$-invariants, by composing $U$ with the group averages $A\circ U:{\mathcal{X}}(\Omega,{\mathcal{C}})\to{\mathcal{C}}^{\prime\prime}$. A natural question is thus whether any ${\mathfrak{G}}$-invariant function can be approximated at arbitrary precision by such a model, for appropriate choices of $B$ and $\sigma$. It is not hard to adapt the standard Universal Approximation Theorems from unstructured vector inputs to show that shallow ‘geometric’ networks are also universal approximators, by properly generalising the group average to a general non-linear invariant. However, as already described in the case of Fourier versus Wavelet invariants, there is a fundamental tension between shallow global invariance and deformation stability. This motivates an alternative representation, which considers instead localised equivariant maps.Assuming that $\Omega$ is further equipped with a distance metric $d$, we call an equivariant map $U$ localised if $(Ux)(u)$ depends only on the values of $x(v)$ for $\mathcal{N}_{u}=\{v:d(u,v)\leq r\}$, for some small radius $r$; the latter set $\mathcal{N}_{u}$ is called the receptive field.

A single layer of local equivariant map $U$ cannot approximate functions with long-range interactions, but a composition of several local equivariant maps $U_{J}\circ U_{J-1}\dots\circ U_{1}$ increases the receptive fieldwhile preserving the stability properties of local equivariants. The receptive field is further increased by interleaving downsampling operators that coarsen the domain (again assuming a metric structure), completing the parallel with Multiresolution Analysis (MRA, see e.g. Mallat (1999)).

Figure 8: Geometric Deep Learning blueprint, exemplified on a graph. A typical Graph Neural Network architecture may contain permutation equivariant layers (computing node-wise features), local pooling (graph coarsening), and a permutation-invariant global pooling layer (readout layer).

In summary, the geometry of the input domain, with knowledge of an underyling symmetry group, provides three key building blocks: (i) a local equivariant map, (ii) a global invariant map, and (iii) a coarsening operator. These building blocks provide a rich function approximation space with prescribed invariance and stability properties by combining them together in a scheme we refer to as the Geometric Deep Learning Blueprint (Figure 8).

Different settings of Geometric Deep Learning

One can make an important distinction between the setting when the domain $\Omega$ is assumed to be fixed and one is only interested in varying input signals defined on that domain, or the domain is part of the input as varies together with signals defined on it. A classical instance of the former case is encountered in computer vision applications, where images are assumed to be defined on a fixed domain (grid). Graph classification is an example of the latter setting, where both the structure of the graph as well as the signal defined on it (e.g. node features) are important. In the case of varying domain, geometric stability (in the sense of insensitivity to the deformation of $\Omega$) plays a crucial role in Geometric Deep Learning architectures.

This blueprint has the right level of generality to be used across a wide range of geometric domains. Different Geometric Deep Learning methods thus differ in their choice of the domain, symmetry group, and the specific implementation details of the aforementioned building blocks. As we will see in the following, a large class of deep learning architectures currently in use fall into this scheme and can thus be derived from common geometric principles.

In the following sections (4.1–4.6) we will describe the various geometric domains focusing on the ‘5G’, and in Sections 5.1–5.8 the specific implementations of Geometric Deep Learning on these domains.

4 Geometric Domains: the 5 Gs

The main focus of our text will be on graphs, grids, groups, geodesics, and gauges. In this context, by ‘groups’ we mean global symmetry transformations in homogeneous space, by ‘geodesics’ metric structures on manifolds, and by ‘gauges’ local reference frames defined on tangent bundles (and vector bundles in general). These notions will be explained in more detail later. In the next sections, we will discuss in detail the main elements in common and the key distinguishing features between these structures and describe the symmetry groups associated with them. Our exposition is not in the order of generality – in fact, grids are particular cases of graphs – but a way to highlight important concepts underlying our Geometric Deep Learning blueprint.

Figure 9: The 5G of Geometric Deep Learning: grids, groups & homogeneous spaces with global symmetry, graphs, geodesics & metrics on manifolds, and gauges (frames for tangent or feature spaces).

4.1 Graphs and Sets

In multiple branches of science, from sociology to particle physics, graphs are used as models of systems of relations and interactions. From our perspective, graphs give rise to a very basic type of invariance modelled by the group of permutations. Furthermore, other objects of interest to us, such as grids and sets, can be obtained as a particular case of graphs.

A graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ is a collection of nodes${\mathcal{V}}$ and edges ${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$ between pairs of nodes. For the purpose of the following discussion, we will further assume the nodes to be endowed with $s$-dimensional node features, denoted by $\mathbf{x}_{u}$ for all $u\in{\mathcal{V}}$. Social networks are perhaps among the most commonly studied examples of graphs, where nodes represent users, edges correspond to friendship relations between them, and node features model user properties such as age, profile picture, etc. It is also often possible to endow the edges, or entire graphs, with features;but as this does not alter the main findings of this section, we will defer discussing it to future work.

The key structural property of graphs is that the nodes in ${\mathcal{V}}$ are usually not assumed to be provided in any particular order, and thus any operations performed on graphs should not depend on the ordering of nodes. The desirable property that functions acting on graphs should satisfy is thus permutation invariance, and it implies that for any two isomorphic graphs, the outcomes of these functions are identical. We can see this as a particular setting of our blueprint, where the domain $\Omega=\mathcal{G}$ and the space $\mathcal{X}(\mathcal{G},\mathbb{R}^{d})$ is that of $d$-dimensional node-wise signals. The symmetry we consider is given by the permutation group $\mathfrak{G}=\Sigma_{n}$, whose elements are all the possible orderings of the set of node indices $\{1,\ldots,n\}$.

Let us first illustrate the concept of permutation invariance on sets, a special case of graphs without edges (i.e., ${\mathcal{E}}=\emptyset$). By stacking the node features as rows of the $n\times d$ matrix $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\top}$, we do effectively specify an ordering of the nodes. The action of the permutation $\mathfrak{g}\in\Sigma_{n}$ on the set of nodes amounts to the reordering of the rows of $\mathbf{X}$, which can be represented as an $n\times n$ permutation matrix $\rho(\mathfrak{g})=\mathbf{P}$,where each row and column contains exactly one $1$ and all the other entries are zeros.

A function $f$ operating on this set is then said to be permutation invariant if, for any such permutation matrix ${\bf P}$, it holds that $f({\bf P}{\bf X})=f({\bf X})$. One simple such function is

$$f({\bf X})=\phi\left(\sum_{u\in{\mathcal{V}}}\psi\left(\mathbf{x}_{u}\right)\right)~{},$$ (6)

where the function $\psi$ is independently applied to every node’s features, and $\phi$ is applied on its sum-aggregated outputs: as sum is independent of the order in which its inputs are provided, such a function is invariant with respect to the permutation of the node set, and is hence guaranteed to always return the same output, no matter how the nodes are permuted.

Functions like the above provide a ‘global’ graph-wise output, but very often, we will be interested in functions that act ‘locally’, in a node-wise manner. For example, we may want to apply some function to update the features in every node, obtaining the set of latent node features. If we stack these latent features into a matrix $\mathbf{H}=\mathbf{F}({\bf X})$is no longer permutation invariant: the order of the rows of ${\bf H}$ should be tied to the order of the rows of ${\bf X}$, so that we know which output node feature corresponds to which input node. We need instead a more fine-grained notion of permutation equivariance, stating that, once we “commit” to a permutation of inputs, it consistently permutes the resulting objects. Formally, $\mathbf{F}(\mathbf{X})$ is a permutation equivariant function if, for any permutation matrix ${\bf P}$, it holds that $\mathbf{F}({\bf P}{\bf X})={\bf P}\mathbf{F}({\bf X})$. A shared node-wise linear transform

$$\mathbf{F}_{\mathbf{\Theta}}({\bf X})={\bf X}{\mathbf{\Theta}}$$ (7)

specified by a weight matrix $\mathbf{\Theta}\in\mathbb{R}^{d\times d^{\prime}}$, is one possible construction of such a permutation equivariant function, producing in our example latent features of the form $\mathbf{h}_{u}=\bm{\Theta}^{\top}\mathbf{x}_{u}$.

This construction arises naturally from our Geometric Deep Learning blueprint. We can first attempt to characterise linear equivariants (functions of the form $\mathbf{F}{\bf PX}=\bf{P}\mathbf{FX}$), for which it is easy to verify that any such map can be written as a linear combination of two generators, the identity $\mathbf{F}_{1}\mathbf{X}={\bf X}$ and the average ${\mathbf{F}_{2}\mathbf{X}}=\frac{1}{n}\bm{1}\bm{1}^{\top}\mathbf{X}=\frac{1}{n}\sum_{u=1}^{n}{\bf x}_{u}$. As will be described in Section 5.4, the popular Deep Sets (Zaheer et al., 2017) architecture follows precisely this blueprint.

We can now generalise the notions of permutation invariance and equivariance from sets to graphs. In the generic setting ${\mathcal{E}}\neq\emptyset$, the graph connectivity can be represented by the $n\times n$ adjacency matrix $\mathbf{A}$,defined as

$$a_{uv}=\begin{cases}1&(u,v)\in{\mathcal{E}}\\ 0&\text{otherwise}.\end{cases}$$ (8)

Note that now the adjacency and feature matrices $\mathbf{A}$ and $\mathbf{X}$ are “synchronised”, in the sense that $a_{uv}$ specifies the adjacency information between the nodes described by the $u$th and $v$th rows of $\mathbf{X}$. Therefore, applying a permutation matrix $\mathbf{P}$ to the node features $\mathbf{X}$ automatically implies applying it to $\mathbf{A}$’s rows and columns, $\mathbf{P}\mathbf{A}\mathbf{P}^{\top}$. We say that (a graph-wise function) $f$ is permutation invariant if

$$f({\bf PX},{\bf PAP}^{\top})=f({\bf X},{\bf A})$$ (9)

and (a node-wise function) $\mathbf{F}$ is permutation equivariant if

$$\mathbf{F}({\bf PX},{\bf PAP}^{\top})={\bf P}\mathbf{F}({\bf X},{\bf A})$$ (10)

for any permutation matrix ${\bf P}$.

Here again, we can first characterise linear equivariant functions.As observed by Maron et al. (2018), any linear $\mathbf{F}$ satisfying equation (10) can be expressed as a linear combination of fifteen linear generators; remarkably, this family of generators is independent of $n$. Amongst these generators, our blueprint specifically advocates for those that are also local, i.e., whereby the output on node $u$ directly depends on its neighbouring nodes in the graph. We can formalise this constraint explicitly in our model construction, by defining what it means for a node to be neighbouring another.

A (undirected) neighbourhood of node $u$, sometimes also called 1-hop, is defined as

$$\mathcal{N}_{u}=\{v:(u,v)\in{\mathcal{E}}\,\mathrm{or}\,(v,u)\in{\mathcal{E}}\}$$ (11)

and the neighbourhood features as the multiset

$$\mathbf{X}_{\mathcal{N}_{u}}=\{\mskip-5.0mu\{\mathbf{x}_{v}:v\in\mathcal{N}_{u}\}\mskip-5.0mu\}.$$ (12)

Operating on 1-hop neighbourhoods aligns well with the locality aspect of our blueprint: namely, defining our metric over graphs as the shortest path distance between nodes using edges in $\mathcal{E}$.

The GDL blueprint thus yields a general recipe for constructing permutation equivariant functions on graphs, by specifying a local function $\phi$ that operates over the features of a node and its neighbourhood, $\phi(\mathbf{x}_{u},\mathbf{X}_{\mathcal{N}_{u}})$. Then, a permutation equivariant function $\mathbf{F}$ can be constructed by applying $\phi$ to every node’s neighbourhood in isolation (see Figure 10):

$$\mathbf{F}({\bf X},{\bf A})=\left[\begin{array}[]{ccc}\rule[2.15277pt]{10.76385pt}{0.5pt}&\phi(\mathbf{x}_{1},\mathbf{X}_{\mathcal{N}_{1}})&\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}&\phi(\mathbf{x}_{2},\mathbf{X}_{\mathcal{N}_{2}})&\rule[2.15277pt]{10.76385pt}{0.5pt}\\ &\vdots&\\ \rule[2.15277pt]{10.76385pt}{0.5pt}&\phi(\mathbf{x}_{n},\mathbf{X}_{\mathcal{N}_{n}})&\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right]$$ (13)

As $\mathbf{F}$ is constructed by applying a shared function $\phi$ to each node locally, its permutation equivariance rests on $\phi$’s output being independent on the ordering of the nodes in $\mathcal{N}_{u}$. Thus, if $\phi$ is built to be permutation invariant, then this property is satisfied. As we will see in future work, the choice of $\phi$ plays a crucial role in the expressive power of such a scheme. When $\phi$ is injective, it is equivalent to one step of the Weisfeiler-Lehman graph isomorphism test, a classical algorithm in graph theory providing a necessary condition for two graphs to be isomorphic by an iterative color refinement procedure.

Figure 10: An illustration of constructing permutation-equivariant functions over graphs, by applying a permutation-invariant function $\phi$ to every neighbourhood. In this case, $\phi$ is applied to the features $\mathbf{x}_{b}$ of node $b$ as well as the multiset of its neighbourhood features, $\mathbf{X}_{\mathcal{N}_{b}}=\{\mskip-5.0mu\{\mathbf{x}_{a},\mathbf{x}_{b},\mathbf{x}_{c},\mathbf{x}_{d},\mathbf{x}_{e}\}\mskip-5.0mu\}$. Applying $\phi$ in this manner to every node’s neighbourhood recovers the rows of the resulting matrix of latents features $\mathbf{H}=\mathbf{F}(\mathbf{X},\mathbf{A})$.

It is also worth noticing that the difference between functions defined on sets and more general graphs in this example is that in the latter case we need to explicitly account for the structure of the domain. As a consequence, graphs stand apart in the sense that the domain becomes part of the input in machine learning problems, whereas when dealing with sets and grids (both particular cases of graphs) we can specify only the features and assume the domain to be fixed. This distinction will be a recurring motif in our discussion. As a result, the notion of geometric stability (invariance to domain deformation) is crucial in most problems of learning on graphs. It straightforwardly follows from our construction that permutation invariant and equivariant functions produce identical outputs on isomorphic (topologically-equivalent) graphs. These results can be generalised to approximately isomorphic graphs, and several results on stability under graph perturbations exist (Levie et al., 2018). We will return to this important point in our discussion on manifolds, which we will use as an vehicle to study such invariance in further detail.

Second, due to their additional structure, graphs and grids, unlike sets, can be coarsened in a non-trivial way, giving rise to a variety of pooling operations.

4.2 Grids and Euclidean spaces

The second type of objects we consider are grids. It is fair to say that the impact of deep learning was particularly dramatic in computer vision, natural language processing, and speech recognition. These applications all share a geometric common denominator: an underlying grid structure. As already mentioned, grids are a particular case of graphs with special adjacency. However, since the order of nodes in a grid is fixed, machine learning models for signals defined on grids are no longer required to account for permutation invariance, and have a stronger geometric prior: translation invariance.

Circulant matrices and Convolutions

Let us dwell on this point in more detail. Assuming for simplicity periodic boundary conditions, we can think of a one-dimensional grid as a ring graphwith nodes indexed by $0,1,\ldots,n-1$ modulo $n$ (which we will omit for notation brevity) and the adjacency matrix with elements $a_{u,u+1\,\mathrm{mod}\,n}=1$ and zero otherwise. There are two main differences from the general graph case we have discussed before. First, each node $u$ has identical connectivity, to its neighbours $u-1$ and $u+1$, and thus structure-wise indistinguishable from the others. Second and more importantly, since the nodes of the grid have a fixed ordering, we also have a fixed ordering of the neighbours: we can call $u-1$ the ‘left neighbour’ and $u+1$ the ‘right neighbour’. If we use our previous recipe for designing a equivariant function $\mathbf{F}$ using a local aggregation function $\phi$, we now have $\mathbf{f}(\mathbf{x}_{u})=\phi(\mathbf{x}_{u-1},\mathbf{x}_{u},\mathbf{x}_{u+1})$ at every node of the grid: $\phi$ does not need to be permutation invariant anymore. For a particular choice of a linear transformation $\phi(\mathbf{x}_{u-1},\mathbf{x}_{u},\mathbf{x}_{u+1})=\theta_{-1}\mathbf{x}_{u-1}+\theta_{0}\mathbf{x}_{u}+\theta_{1}\mathbf{x}_{u+1}$, we can write $\mathbf{F}(\mathbf{X})$ as a matrix product,

$$\mathbf{F}(\mathbf{X})=\left[\begin{array}[]{ccccc}\theta_{0}&\theta_{1}&&&\theta_{-1}\\ \theta_{-1}&\theta_{0}&\theta_{1}&&\\ &\ddots&\ddots&\ddots&\\ &&\theta_{-1}&\theta_{0}&\theta_{1}\\ \theta_{1}&&&\theta_{-1}&\theta_{0}\par\end{array}\right]\left[\begin{array}[]{ccc}\rule[2.15277pt]{10.76385pt}{0.5pt}&\mathbf{x}_{0}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}&\mathbf{x}_{1}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\ &\vdots&\\ \rule[2.15277pt]{10.76385pt}{0.5pt}&\mathbf{x}_{n-2}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}&\mathbf{x}_{n-1}&\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right]$$

Note this very special multi-diagonal structure with one element repeated along each diagonal, sometimes referred to as “weight sharing” in the machine learning literature.

More generally, given a vector $\bm{\theta}=(\theta_{0},\ldots,\theta_{n-1})$, a circulant matrix $\mathbf{C}(\bm{\theta})=(\theta_{u-v\,\mathrm{mod}\,n})$ is obtained by appending circularly shifted versions of the vector $\bm{\theta}$. Circulant matrices are synonymous with discrete convolutions,

$$(\mathbf{x}\star\bm{\theta})_{u}=\sum_{v=0}^{n-1}x_{v\,\mathrm{mod}\,n}\,\,\theta_{u-v\,\mathrm{mod}\,n}$$

as one has $\mathbf{C}(\bm{\theta})\mathbf{x}=\mathbf{x}\star\bm{\theta}$. A particular choice of $\bm{\theta}=(0,1,0,\ldots,0)^{\top}$ yields a special circulant matrix that shifts vectors to the right by one position. This matrix is called the (right) shift or translation operator and denoted by $\mathbf{S}$.

Circulant matrices can be characterised by their commutativity property: the product of circulant matrices is commutative, i.e. $\mathbf{C}(\bm{\theta})\mathbf{C}(\bm{\eta})=\mathbf{C}(\bm{\eta})\mathbf{C}(\bm{\theta})$ for any $\bm{\theta}$ and $\bm{\eta}$. Since the shift is a circulant matrix, we get the familiar translation or shift equivariance of the convolution operator,

$$\mathbf{S}\mathbf{C}(\bm{\theta})\mathbf{x}=\mathbf{C}(\bm{\theta})\mathbf{S}\mathbf{x}.$$

Such commutativity property should not be surprising, since the underlying symmetry group (the translation group) is Abelian. Moreover, the opposite direction appears to be true as well, i.e. a matrix is circulant iff it commutes with shift. This, in turn, allows us to define convolution as a translation equivariant linear operation, and is a nice illustration of the power of geometric priors and the overall philosophy of Geometric ML: convolution emerges from the first principle of translational symmetry.

Note that unlike the situation on sets and graphs, the number of linearly independent shift-equivariant functions (convolutions) grows with the size of the domain (since we have one degree of freedom in each diagonal of a circulant matrix). However, the scale separation prior guarantees filters can be local, resulting in the same $\Theta(1)$-parameter complexity per layer, as we will verify in Section 5.1 when discussing the use of these principles in the implementation of Convolutional Neural Network architectures.

Derivation of the discrete Fourier transform

We have already mentioned the Fourier transform and its connection to convolution: the fact that the Fourier transform diagonalises the convolution operation is an important property used in signal processing to perform convolution in the frequency domain as an element-wise product of the Fourier transforms. However, textbooks usually only state this fact, rarely explaining where the Fourier transform comes from and what is so special about the Fourier basis. Here we can show it, demonstrating once more how foundational are the basic principles of symmetry.

For this purpose, recall a fact from linearalgebra that (diagonalisable) matrices are joinly diagonalisable iff they mutually commute. In other words, there exists a common eigenbasis for all the circulant matrices, in which they differ only by their eigenvalues. We can therefore pick one circulant matrix and compute its eigenvectors—we are assured that these will be the eigenvectors of all other circulant matrices as well. It is convenient to pick the shift operator, for which the eigenvectors happen to be the discrete Fourier basis

$$\bm{\varphi}_{k}=\frac{1}{\sqrt{n}}\left(1,e^{\frac{2\pi\mathrm{i}k}{n}},e^{\frac{4\pi\mathrm{i}k}{n}},\ldots,e^{\frac{2\pi\mathrm{i}(n-1)k}{n}}\right)^{\top},\hskip 14.22636ptk=0,1,\ldots,n-1,$$

which we can arrange into an $n\times n$ Fourier matrix $\bm{\Phi}=(\bm{\varphi}_{0},\ldots,\bm{\varphi}_{n-1})$. Multiplication by $\bm{\Phi}^{*}$gives the Discrete Fourier Transform (DFT), and by $\bm{\Phi}$ the inverse DFT,

$$\hat{x}_{k}=\frac{1}{\sqrt{n}}\sum_{u=0}^{n-1}x_{u}e^{-\frac{2\pi\mathrm{i}ku}{n}}\hskip 42.67912pt{x}_{u}=\frac{1}{\sqrt{n}}\sum_{k=0}^{n-1}\hat{x}_{k}e^{+\frac{2\pi\mathrm{i}ku}{n}}.$$

Since all circulant matrices are jointly diagonalisable,they are also diagonalised by the Fourier transform and differ only in their eigenvalues. Since the eigenvalues of the circulant matrix $\mathbf{C}(\bm{\theta})$ are the Fourier transform of the filter (see e.g. Bamieh (2018)), $\hat{\bm{\theta}}=\bm{\Phi}^{*}\bm{\theta}$, we obtain the Convolution Theorem:

$$\mathbf{C}(\bm{\theta})\mathbf{x}=\bm{\Phi}\left[\begin{array}[]{ccc}\hat{\theta}_{0}&&\\ &\ddots&\\ &&\hat{\theta}_{n-1}\end{array}\right]\bm{\Phi}^{*}\mathbf{x}=\bm{\Phi}(\hat{\bm{\theta}}\odot\hat{\mathbf{x}})$$

Because the Fourier matrix $\bm{\Phi}$ has a special algebraic structure, the products $\bm{\Phi}^{\star}\mathbf{x}$ and $\bm{\Phi}\mathbf{x}$ can be computed with $\mathcal{O}(n\log n)$ complexity using a Fast Fourier Transform (FFT) algorithm. This is one of the reasons why frequency-domain filtering is so popular in signal processing; furthermore, the filter is typically designed directly in the frequency domain, so the Fourier transform $\hat{\bm{\theta}}$ is never explicitly computed.

Besides the didactic value of the derivation of the Fourier transform and convolution we have done here, it provides a scheme to generalise these concepts to graphs. Realising that the adjacency matrix of the ring graph is exactly the shift operator, one can can develop the graph Fourier transform and an analogy of the convolution operator by computing the eigenvectors of the adjacency matrix (see e.g. Sandryhaila and Moura (2013)). Early attempts to develop graph neural networks by analogy to CNNs, sometimes termed ‘spectral GNNs’, exploited this exact blueprint.We will see in Sections 4.4–4.6 that this analogy has some important limitations. The first limitation comes from the fact that a grid is fixed, and hence all signals on it can be represented in the same Fourier basis. In contrast, on general graphs, the Fourier basis depends on the structure of the graph. Hence, we cannot directly compare Fourier transforms on two different graphs — a problem that translated into a lack of generalisation in machine learning problems. Secondly, multi-dimensional grids, which are constructed as tensor products of one-dimensional grids, retain the underlying structure: the Fourier basis elements and the corresponding frequencies (eigenvalues) can be organised in multiple dimensions. In images, for example, we can naturally talk about horizontal and vertical frequency and filters have a notion of direction. On graphs, the structure of the Fourier domain is one-dimensional, as we can only organise the Fourier basis functions by the magnitude of the corresponding frequencies. As a result, graph filters are oblivious of direction or isotropic.

Derivation of the continuous Fourier transform

For the sake of completeness, and as a segway for the next discussion, we repeat our analysis in the continuous setting. Like in Section 3.4, consider functions defined on $\Omega=\mathbb{R}$ and the translation operator $(S_{v}f)(u)=f(u-v)$ shifting $f$ by some position $v$. Applying $S_{v}$ to the Fourier basis functions $\varphi_{\xi}(u)=e^{\mathrm{i}\xi u}$ yields, by associativity of the exponent,

$$S_{v}e^{\mathrm{i}\xi u}=e^{\mathrm{i}\xi(u-v)}=e^{-\mathrm{i}\xi v}e^{\mathrm{i}\xi u},$$

i.e., $\varphi{u}_{\xi}(u)$ is the complex eigenvector of $S_{v}$ with the complex eigenvalue $e^{-\mathrm{i}\xi v}$ – exactly mirroring the situation we had in the discrete setting. Since $S_{v}$ is a unitary operator (i.e., $\|S_{v}x\|_{p}=\|x\|_{p}$ for any $p$ and $x\in L_{p}(\mathbb{R})$), any eigenvalue $\lambda$ must satisfy $|\lambda|=1$, which corresponds precisely to the eigenvalues $e^{-i\xi v}$ found above. Moreover, the spectrum of the translation operator is simple, meaning that two functions sharing the same eigenvalue must necessarily be collinear. Indeed, suppose that $S_{v}f=e^{-\mathrm{i}\xi_{0}v}f$ for some $\xi_{0}$. Taking the Fourier transform in both sides, we obtain

$$\forall~{}\xi~{},~{}e^{-\mathrm{i}\xi v}\hat{f}(\xi)=e^{-\mathrm{i}\xi_{0}v}\hat{f}(\xi)~{},$$

which implies that $\hat{f}(\xi)=0$ for $\xi\neq\xi_{0}$, thus $f=\alpha\varphi_{\xi_{0}}$.

For a general linear operator $C$ that is translation equivariant ($S_{v}C=CS_{v}$), we have

$$S_{v}Ce^{\mathrm{i}\xi u}=CS_{v}e^{\mathrm{i}\xi u}=e^{-\mathrm{i}\xi v}Ce^{\mathrm{i}\xi u},$$

implying that $Ce^{\mathrm{i}\xi u}$ is also an eigenfunctionof $S_{v}$ with eigenvalue $e^{-\mathrm{i}\xi v}$, from where it follows from the simplicity of spectrum that $Ce^{\mathrm{i}\xi u}=\beta\varphi_{\xi}(u)$; in other words, the Fourier basis is the eigenbasis of all translation equivariant operators. As a result, $C$ is diagonal in the Fourier domain and can be expressed as $Ce^{\mathrm{i}\xi u}=\hat{p}_{C}(\xi)e^{\mathrm{i}\xi u}$, where $\hat{p}_{C}(\xi)$ is a transfer function acting on different frequencies $\xi$. Finally, for an arbitrary function $x(u)$, by linearity,

$$(Cx)(u) = C\int_{-\infty}^{+\infty}\hat{x}(\xi)e^{\mathrm{i}\xi u}\mathrm{d}\xi=\int_{-\infty}^{+\infty}\hat{x}(\xi)\hat{p}_{C}(\xi)e^{\mathrm{i}\xi u}\mathrm{d}\xi$$

$$= \int_{-\infty}^{+\infty}p_{C}(v)x(u-v)\mathrm{d}v~{}=(x\star p_{C})(u),$$

where $p_{C}(u)$ is the inverse Fourier transform of $\hat{p}_{C}(\xi)$. It thus follows that every linear translation equivariant operator is a convolution.

4.3 Groups and Homogeneous spaces

Our discussion of grids highlighted how shifts and convolutions are intimately connected: convolutions are linear shift-equivariantoperations, and vice versa, any shift-equivariant linear operator is a convolution. Furthermore, shift operators can be jointly diagonalised by the Fourier transform. As it turns out, this is part of a far larger story: both convolution and the Fourier transform can be defined for any group of symmetries that we can sum or integrate over.

Consider the Euclidean domain $\Omega=\mathbb{R}$. We can understand the convolution as a pattern matching operation: we match shifted copies of a filter $\theta(u)$ with an input signal $x(u)$. The value of the convolution $(x\star\theta)(u)$ at a point $u$ is the inner product of the signal $x$ with the filter shifted by $u$,

$$(x\star\theta)(u)=\langle x,S_{u}\theta\rangle=\int_{\mathbb{R}}x(v)\theta(u+v)\mathrm{d}v.$$

Note that in this case $u$ is both a point on the domain $\Omega=\mathbb{R}$ and also an element of the translation group, which we can identify with the domain itself, ${\mathfrak{G}}=\mathbb{R}$. We will now show how to generalise this construction, by simply replacing the translation group by another group ${\mathfrak{G}}$ acting on $\Omega$.

Group convolution

As discussed in Section 3, the action of the group ${\mathfrak{G}}$ on the domain $\Omega$ induces a representation $\rho$ of ${\mathfrak{G}}$ on the space of signals $\mathcal{X}(\Omega)$ via $\rho({\mathfrak{g}})x(u)=x({\mathfrak{g}}^{-1}u)$. In the above example, ${\mathfrak{G}}$ is the translation group whose elements act by shifting the coordinates, $u+v$, whereas $\rho({\mathfrak{g}})$ is the shift operator acting on signals as $(S_{v}x)(u)=x(u-v)$. Finally, in order to apply a filter to the signal, we invoke our assumption of $\mathcal{X}(\Omega)$ being a Hilbert space, with an inner product

$$\langle x,\theta\rangle=\int_{\Omega}x(u)\theta(u)\mathrm{d}u,$$

where we assumed, for the sake of simplicity, scalar-valued signals, $\mathcal{X}(\Omega,\mathbb{R})$; in general the inner product has the form of equation (2).

Having thus defined how to transform signals and match them with filters, we can define the group convolution for signals on $\Omega$,

$$(x\star\theta)({\mathfrak{g}})=\langle x,\rho({\mathfrak{g}})\theta\rangle=\int_{\Omega}x(u)\theta({\mathfrak{g}}^{-1}u)\mathrm{d}u.$$ (14)

Note that $x\star\theta$ takes values on the elements ${\mathfrak{g}}$ of our group ${\mathfrak{G}}$ rather than points on the domain $\Omega$. Hence, the next layer, which takes $x\star\theta$ as input, should act on signals defined on to the group ${\mathfrak{G}}$, a point we will return to shortly.

Just like how the traditional Euclidean convolution is shift-equivariant, the more general group convolution is ${\mathfrak{G}}$-equivariant. The key observation is that matching the signal $x$ with a ${\mathfrak{g}}$-transformed filter $\rho({\mathfrak{g}})\theta$ is the same as matching the inverse transformed signal $\rho({\mathfrak{g}}^{-1})x$ with the untransformed filter $\theta$. Mathematically, this can be expressed as $\langle x,\rho({\mathfrak{g}})\theta\rangle=\langle\rho({\mathfrak{g}}^{-1})x,\theta\rangle$. With this insight, ${\mathfrak{G}}$-equivariance of the group convolution (14) follows immediately from its definition and the defining property $\rho({\mathfrak{h}}^{-1})\rho({\mathfrak{g}})=\rho({\mathfrak{h}}^{-1}{\mathfrak{g}})$ of group representations,

$$(\rho({\mathfrak{h}})x\star\theta)({\mathfrak{g}})=\langle\rho({\mathfrak{h}})x,\rho({\mathfrak{g}})\theta\rangle=\langle x,\rho({\mathfrak{h}}^{-1}{\mathfrak{g}})\theta\rangle=\rho({\mathfrak{h}})(x\star\theta)({\mathfrak{g}}).$$

Let us look at some examples. The case of one-dimensional grid we have studied above is obtained with the choice $\Omega=\mathbb{Z}_{n}=\{0,\ldots,n-1\}$ and the cyclic shift group ${\mathfrak{G}}=\mathbb{Z}_{n}$. The group elements in this case are cyclic shifts of indices, i.e., an element ${\mathfrak{g}}\in{\mathfrak{G}}$ can be identified with some $u=0,\ldots,n-1$ such that ${\mathfrak{g}}.v=v-u\,\mathrm{mod}\,n$, whereas the inverse element is ${\mathfrak{g}}^{-1}.v=v+u\,\mathrm{mod}\,n$. Importantly, in this example the elements of the group (shifts) are also elements of the domain (indices). We thus can, with some abuse of notation, identify the two structures (i.e., $\Omega={\mathfrak{G}}$); our expression for the group convolution in this case

$$(x\star\theta)({\mathfrak{g}})=\sum_{v=0}^{n-1}x_{v}\,\theta_{{\mathfrak{g}}^{-1}v},$$

leads to the familiar convolution $\displaystyle(x\star\theta)_{u}=\sum_{v=0}^{n-1}x_{v}\,\theta_{v+u\,\,\mathrm{mod}\,n}$.

Spherical convolution

Now consider the two-dimensional sphere $\Omega=\mathbb{S}^{2}$ with the group of rotations, the special orthogonal group ${\mathfrak{G}}=\mathrm{SO}(3)$. While chosen for pedagogical reason, this example is actually very practical and arises in numerous applications. In astrophysics, for example, observational data often naturally has spherical geometry. Furthermore, spherical symmetries are very important in applications in chemistry when modeling molecules and trying to predict their properties, e.g. for the purpose of virtual drug screening.

Representing a point on the sphere as a three-dimensional unit vector $\mathbf{u}:\|\mathbf{u}\|=1$, the action of the group can be represented as a $3\times 3$ orthogonal matrix $\mathbf{R}$ with $\mathrm{det}(\mathbf{R})=1$. The spherical convolution can thus be written as the inner product between the signal and the rotated filter,

$$(x\star\theta)(\mathbf{R})=\int_{\mathbb{S}^{2}}x(\textbf{u})\theta(\mathbf{R}^{-1}\mathbf{u})\mathrm{d}\mathbf{u}.$$

The first thing to note is than now the group is not identical to the domain: the group $\mathrm{SO}(3)$ is a Lie group that is in fact a three-dimensional manifold, whereas $\mathbb{S}^{2}$ is a two-dimensional one. Consequently, in this case, unlike the previous example, the convolution is a function on $\operatorname{SO}(3)$ rather than on $\Omega$.

This has important practical consequences: in our Geometric Deep Learning blueprint, we concatenate multiple equivariant maps (“layers” in deep learning jargon) by applying a subsequent operator to the output of the previous one. In the case of translations, we can apply multiple convolutions in sequence, since their outputs are all defined on the same domain $\Omega$. In the general setting, since $x\star\theta$ is a function on ${\mathfrak{G}}$ rather than on $\Omega$, we cannot use exactly the same operation subsequently—it means that the next operation has to deal with signals on ${\mathfrak{G}}$, i.e. $x\in\mathcal{X}({\mathfrak{G}})$. Our definition of group convolution allows this case: we take as domain $\Omega={\mathfrak{G}}$ acted on by ${\mathfrak{G}}$ itself via the group action $({\mathfrak{g}},{\mathfrak{h}})\mapsto{\mathfrak{g}}{\mathfrak{h}}$ defined by the composition operation of ${\mathfrak{G}}$. This yields the representation $\rho({\mathfrak{g}})$ acting on $x\in\mathcal{X}({\mathfrak{G}})$ by $(\rho({\mathfrak{g}})x)({\mathfrak{h}})=x({\mathfrak{g}}^{-1}{\mathfrak{h}})$. Just like before, the inner product is defined by integrating the point-wise product of the signal and the filter over the domain, which now equals $\Omega={\mathfrak{G}}$. In our example of spherical convolution, a second layer of convolution would thus have the form

$$((x\star\theta)\star\phi)(\mathbf{R})=\int_{\mathrm{SO}(3)}(x\star\theta)(\textbf{Q})\phi(\mathbf{R}^{-1}\mathbf{Q})\mathrm{d}\mathbf{Q}.$$

Since convolution involves inner product that in turn requires integrating over the domain $\Omega$, we can only use it on domains $\Omega$ that are small (in the discrete case) or low-dimensional (in the continuous case). For instance, we can use convolutions on the plane $\mathbb{R}^{2}$ (two dimensional) or special orthogonal group $\operatorname{SE}(3)$ (three dimensional), or on the finite set of nodes of a graph ($n$-dimensional), but we cannot in practice perform convolution on the group of permutations $\Sigma_{n}$, which has $n!$ elements. Likewise, integrating over higher-dimensional groups like the affine group (containing translations, rotations, shearing and scaling, for a total of $6$ dimensions) is not feasible in practice. Nevertheless, as we have seen in Section 5.3, we can still build equivariant convolutions for large groups ${\mathfrak{G}}$ by working with signals defined on low-dimensional spaces $\Omega$ on which ${\mathfrak{G}}$ acts. Indeed, it is possible to show that any equivariant linear map $f:\mathcal{X}(\Omega)\rightarrow\mathcal{X}(\Omega^{\prime})$ between two domains $\Omega,\Omega^{\prime}$ can be written as a generalised convolution similar to the group convolution discussed here.

Second, we note that the Fourier transform we derived in the previous section from the shift-equivariance property of the convolution can also be extended to a more general case by projecting the signal onto the matrix elements of irreducible representations of the symmetry group. We will discuss this in future work. In the case of $\operatorname{SO}(3)$ studied here, this gives rise to the spherical harmonics and Wigner D-functions, which find wide applications in quantum mechanics and chemistry.

Finally, we point to the assumption that has so far underpinned our discussion in this section: whether $\Omega$ was a grid, plane, or the sphere, we could transform every point into any other point, intuitively meaning that all the points on the domain “look the same.” A domain $\Omega$ with such property is called a homogeneous space, where for any $u,v\in\Omega$ there exists $\mathfrak{g}\in{\mathfrak{G}}$ such that $\mathfrak{g}.u=v$. In the next section we will try to relax this assumption.

4.4 Geodesics and Manifolds

In our last example, the sphere $\mathbb{S}^{2}$ was a manifold, albeit a special one with a global symmetry group due to its homogeneous structure. Unfortunately, this is not the case for the majority of manifolds, which typically do not have global symmetries. In this case, we cannot straightforwardly define an action of ${\mathfrak{G}}$ on the space of signals on $\Omega$ and use it to ‘slide’ filters around in order to define a convolution as a direct generalisation of the classical construction. Nevertheless, manifolds do have two types of invariance that we will explore in this section: transformations preserving metric structure and local reference frame change.

While for many machine learning readers manifolds might appear as somewhat exotic objects, they are in fact very common in various scientific domains. In physics, manifolds play a central role as the model of our Universe — according to Einstein’s General Relativity Theory, gravity arises from the curvature of the space-time, modeled as a pseudo-Riemannian manifold. In more ‘prosaic’ fields such as computer graphics and vision, manifolds are a common mathematical model of 3D shapes.The broad spectrum of applications of such models ranges from virtual and augmented reality and special effects obtained by means of ‘motion capture’ to structural biology dealing with protein interactions that stick together (‘bind’ in chemical jargon) like pieces of 3D puzzle. The common denominator of these applications is the use of a manifold to represent the boundary surface of some 3D object.

There are several reasons why such models are convenient.First, they offer a compact description of the 3D object, eliminating the need to allocate memory to ‘empty space’ as is required in grid-based representations. Second, they allow to ignore the internal structure of the object. This is a handy property for example in structural biology where the internal folding of a protein molecule is often irrelevant for interactions that happen on the molecular surface. Third and most importantly, one often needs to deal with deformable objects that undergo non-rigid deformations. Our own body is one such example, and many applications in computer graphics and vision, such as the aforementioned motion capture and virtual avatars, require deformation invariance. Such deformations can be modelled very well as transformations that preserve the intrinsic structure of a (Riemannian) manifold, namely the distances between points measured along the manifold, without regard to the way the manifold is embedded in the ambient space.

We should emphasise that manifolds fall under the setting of varying domains in our Geometric Deep Learning blueprint, and in this sense are similar to graphs. We will highlight the importance of the notion of invariance to domain deformations – what we called ‘geometric stability’ in Section 3.3. Since differential geometry is perhaps less familiar to the machine learning audience, we will introduce the basic concepts required for our discussion and refer the reader to Penrose (2005) for their detailed exposition.

Riemannian manifolds

Since the formal definition of a manifoldis somewhat involved, we prefer to provide an intuitive picture at the expense of some precision. In this context, we can think of a (differentiable or smooth) manifold as a smooth multidimensional curved surface that is locally Euclidean, in the sense that any small neighbourhood around any point it can be deformed to a neighbourhood of $\mathbb{R}^{s}$; in this case the manifold is said to be $s$-dimensional. This allows us to locally approximate the manifold around point $u$ through the tangent space $T_{u}\Omega$. The latter can be visualised by thinking of a prototypical two-dimensional manifold, the sphere, and attaching a plane to it at a point: with sufficient zoom, the spherical surface will seem planar (Figure 11). The collection of all tangent spaces is called the tangent bundle, denoted $T\Omega$; we will dwell on the concept of bundles in more detail in Section 4.5.

A tangent vector, which we denote by $X\in T_{u}\Omega$, can be thought of as a local displacement from point $u$. In order to measure the lengths of tangent vectors and angles between them, we need to equip the tangent space with additional structure, expressed as a positive-definite bilinear function $g_{u}:T_{u}\Omega\times T_{u}\Omega\rightarrow\mathbb{R}$ depending smoothly on $u$. Such a function is called a Riemannian metric, in honour of Bernhardt Riemann who introduced the concept in 1856, and can be thought of as an inner product on the tangent space, $\langle X,Y\rangle_{u}=g_{u}(X,Y)$, which is an expression of the angle between any two tangent vectors $X,Y\in T_{u}\Omega$. The metric also induces a norm $\|X\|_{u}=g_{u}^{1/2}(X,X)$ allowing to locally measure lengths of vectors.

We must stress that tangent vectors are abstract geometric entities that exists in their own right and are coordinate-free. If we are to express a tangent vector $X$ numerically as an array of numbers, we can only represent it as a list of coordinates $\mathbf{x}=(x_{1},\ldots,x_{s})$ relative to some local basis$\{X_{1},\ldots X_{s}\}\subseteq T_{u}\Omega$. Similarly, the metric can be expressed as an $s\times s$ matrix $\mathbf{G}$ with elements $g_{ij}=g_{u}(X_{i},X_{j})$ in that basis. We will return to this point in Section 4.5.

Figure 11: Basic notions of Riemannian geometry illustrated on the example of the two-dimensional sphere $\mathbb{S}^{2}=\{\mathbf{u}\in\mathbb{R}^{3}:\|\mathbf{u}\|=1\}$, realised a subset (sub-manifold) of $\mathbb{R}^{3}$. The tangent space to the sphere is given as $T_{\mathbf{u}}\mathbb{S}^{2}=\{\mathbf{x}\in\mathbb{R}^{3}:\mathbf{x}^{\top}\mathbf{u}=0\}$ and is a 2D plane – hence this is a 2-dimensional manifold. The Riemannian metric is simply the Euclidean inner product restricted to the tangent plane, $\langle\mathbf{x},\mathbf{y}\rangle_{\mathbf{u}}=\mathbf{x}^{\top}\mathbf{y}$ for any $\mathbf{x},\mathbf{x}\in T_{\mathbf{u}}\mathbb{S}^{2}$. The exponential map is given by $\exp_{\mathbf{u}}(\mathbf{x})=\cos(\|\mathbf{x}\|)\mathbf{u}+\frac{\sin(\|\mathbf{x}\|)}{\|\mathbf{x}\|}\mathbf{x}$, for $\mathbf{x}\in T_{\mathbf{u}}\mathbb{S}^{2}$. Geodesics are great arcs of length $d(\mathbf{u},\mathbf{v})=\cos^{-1}(\mathbf{u}^{\top}\mathbf{v})$.

A manifold equipped with a metric is called a Riemannian manifold and properties that can be expressed entirely in terms of the metric are said to be intrinsic. This is a crucial notion for our discussion, as according to our template, we will be seeking to construct functions acting on signals defined on $\Omega$ that are invariant to metric-preserving transformations called isometries that deform the manifold without affecting its local structure.If such functions can be expressed in terms of intrinsic quantities, they are automatically guaranteed to be isometry-invariant and thus unaffected by isometric deformations. These results can be further extended to dealing with approximate isometries; this is thus an instance of the geometric stability (domain deformation) discussed in our blueprint.

While, as we noted, the definition of a Riemannian manifold does not require a geometric realisation in any space, it turns out that any smooth Riemannian manifold can be realised as a subset of a Euclidean space of sufficiently high dimension (in which case it is said to be ‘embedded’ in that space) by using the structure of the Euclidean space to induce a Riemannian metric. Such an embedding is however not necessarily unique – as we will see, two different isometric realisations of a Riemannian metric are possible.

Scalar and Vector fields

Since we are interested in signals defined on $\Omega$, we need to provide the proper notion of scalar- and vector-valued functions on manifolds. A (smooth) scalar field is a function of the form $x:\Omega\rightarrow\mathbb{R}$. Scalar fields form a vector space $\mathcal{X}(\Omega,\mathbb{R})$ that can be equipped with the inner product

$$\langle x,y\rangle=\int_{\Omega}x(u)y(u)\mathrm{d}u,$$ (15)

where $\mathrm{d}u$ is the volume element induced by the Riemannian metric. A (smooth) tangent vector field is a function of the form $X:\Omega\rightarrow T\Omega$ assigning to each point a tangent vector in the respective tangent space, $u\mapsto X(u)\in T_{u}\Omega$. Vector fieldsalso form a vector space $\mathcal{X}(\Omega,T\Omega)$ with the inner product defined through the Riemannian metric,

$$\langle X,Y\rangle=\int_{\Omega}g_{u}(X(u),Y(u))\mathrm{d}u.$$ (16)

Intrinsic gradient

Another way to think of (and actually define) vector fields is as a generalised notion of derivative. In classical calculus, one can locally linearise a (smooth) function through the differential $\mathrm{d}x(u)=x(u+\mathrm{d}u)-x(u)$, which provides the change of the value of the function $x$ at point $u$ as a result of an inifinitesimal displacement $\mathrm{d}u$. However, in our case the naïve use of this definition is impossible, since expressions of the form “$u+\mathrm{d}u$” are meaningless on manifolds due to the lack of a global vector space structure.

The solution is to use tangent vectors as a model of local infinitesimal displacement. Given a smooth scalar field $x\in\mathcal{X}(\Omega,\mathbb{R})$, we can think of a (smooth) vector field as a linear map $Y:\mathcal{X}(\Omega,\mathbb{R})\rightarrow\mathcal{X}(\Omega,\mathbb{R})$ satisfying the properties of a derivation: $Y(c)=0$ for any constant $c$ (corresponding to the intuition that constant functions have vanishing derivatives), $Y(x+z)=Y(x)+Y(z)$ (linearity), and $Y(xz)=Y(x)z+xY(z)$ (product or Leibniz rule), for any smooth scalar fields $x,z\in\mathcal{X}(\Omega,\mathbb{R})$. It can be shown that one can use these properties to define vector fields axiomatically. The differential $\mathrm{d}x(Y)=Y(x)$ can be viewed as an operator $(u,Y)\mapsto Y(x)$ and interpreted as follows: the change of $x$ as the result of displacement $Y\in T_{u}\Omega$ at point $u$ is given by $\mathrm{d}_{u}x(Y)$. It is thus an extension of the classical notion of directional derivative.

Alternatively, at each point $u$ the differential can be regarded as a linear functional $\mathrm{d}x_{u}:T_{u}\Omega\rightarrow\mathbb{R}$ acting on tangent vectors $X\in T_{u}\Omega$. Linear functionals on a vector space are called dual vectors or covectors; if in addition we are given an inner product (Riemannian metric), a dual vector can always be represented as

$$\mathrm{d}x_{u}(X)=g_{u}(\nabla x(u),X).$$

The representation of the differential at point $u$ is a tangent vector $\nabla x(u)\in T_{u}\Omega$ called the (intrinsic) gradient of $x$; similarly to the gradient in classical calculus, it can be thought of as the direction of the steepest increase of $x$. The gradient considered as an operator $\nabla:\mathcal{X}(\Omega,\mathbb{R})\rightarrow\mathcal{X}(\Omega,T\Omega)$ assigns at each point $x(u)\mapsto\nabla x(u)\in T_{u}\Omega$; thus, the gradient of a scalar field $x$ is a vector field $\nabla x$.

Geodesics

Now consider a smooth curve $\gamma:[0,T]\rightarrow\Omega$ on the manifold with endpoints $u=\gamma(0)$ and $v=\gamma(T)$. The derivative of the curve at point $t$ is a tangent vector $\gamma^{\prime}(t)\in T_{\gamma(t)}\Omega$ called the velocity vector. Among all the curves connecting points $u$ and $v$, we are interested in those of minimum length, i.e., we are seeking $\gamma$ minimising the length functional

$$\ell(\gamma)=\int_{0}^{T}\|\gamma^{\prime}(t)\|_{\gamma(t)}\mathrm{d}t=\int_{0}^{T}g_{\gamma(t)}^{1/2}(\gamma^{\prime}(t),\gamma^{\prime}(t))\mathrm{d}t.$$

Such curves are called geodesics (from the Greek γεοδαιςία, literally ‘division of Earth’) and they play important role in differential geometry. Crucially to our discussion, the way we defined geodesics is intrinsic, as they depend solely on the Riemannian metric (through the length functional).

Readers familiar with differential geometry might recall that geodesics are a more general concept and their definition in fact does not necessarily require a Riemannian metric but a connection (also called a covariant derivative, as it generalises the notion of derivative to vector and tensor fields), which is defined axiomatically, similarly to our construction of the differential. Given a Riemannian metric, there exists a unique special connection called theLevi-Civita connection which is often tacitly assumed in Riemannian geometry. Geodesics arising from this connection are the length-minimising curves we have defined above.

We will show next how to use geodesics to define a way to transport tangent vectors on the manifold (parallel transport), create local intrinsic maps from the manifold to the tangent space (exponential map), and define distances (geodesic metric). This will allow us to construct convolution-like operations by applying a filter locally in the tangent space.

Parallel transport

One issue we have already encountered when dealing with manifolds is that we cannot directly add or subtract two points $u,v\in\Omega$. The same problem arises when trying to compare tangent vectors at different points: though they have the same dimension, they belong to different spaces, e.g. $X\in T_{u}\Omega$ and $Y\in T_{v}\Omega$, and thus not directly comparable. Geodesics provide a mechanism to move vectors from one point to another, in the following way: let $\gamma$ be a geodesic connecting points $u=\gamma(0)$ and $v=\gamma(T)$ and let $X\in T_{u}\Omega$. We can define a new set of tangent vectors along the geodesic, $X(t)\in T_{\gamma(t)}\Omega$ such that the length of $X(t)$ and the angle (expressed through the Riemannian metric) between it and the velocity vector of the curve is constant,

$$g_{\gamma(t)}(X(t),\gamma^{\prime}(t))=g_{u}(X,\gamma^{\prime}(0))=\mathrm{const},\quad\quad\|X(t)\|_{\gamma(t)}=\|X\|_{u}=\mathrm{const}.$$

As a result, we get a unique vector $X(T)\in T_{v}\Omega$ at the end point $v$.

The map $\Gamma_{u\rightarrow v}(X):T_{u}\Omega\rightarrow T_{u}\Omega$ and $T_{v}\Omega$ defined as $\Gamma_{u\rightarrow v}(X)=X(T)$ using the above notation is called parallel transport or connection; the latter term implying it is a mechanism to ‘connect’ between the tangent spaces $T_{u}\Omega$ and $T_{v}\Omega$. Due to the angle and length preservation conditions, parallel transport amounts to only rotation of the vector, so it can be associated with an element of the special orthogonal group $\mathrm{SO}(s)$ (called the structure group of the tangent bundle),which we will denote by ${\mathfrak{g}}_{u\rightarrow v}$ and discuss in further detail in Section 4.5.

As we mentioned before, a connection can be defined axiomatically independently of the Riemannian metric, providing thus an abstract notion of parallel transport along any smooth curve. The result of such transport, however, depends on the path taken.

Exponential map

Locally around a point $u$, it is always possible to define a unique geodesic in a given direction $X\in T_{u}\Omega$, i.e. such that $\gamma(0)=u$ and $\gamma^{\prime}(0)=X$. When $\gamma_{X}(t)$ is defined for all $t\geq 0$ (that is, we can shoot the geodesic from a point $u$ for as long as we like), the manifold is said to be geodesically complete and the exponential map is defined on the whole tangent space. Since compact manifolds are geodesically complete, we can tacitly assume this convenient property.

This definition of geodesic provided a point and a direction gives a natural mapping from (a subset of) the tangent space $T_{u}\Omega$ to $\Omega$ called the exponential map$\exp:B_{r}(0)\subset T_{u}\Omega\rightarrow\Omega$, which is defined by taking a unit step along the geodesic in the direction $X$, i.e., $\exp_{u}(X)=\gamma_{X}(1)$. The exponential map $\exp_{u}$ is a local diffeomorphism, as it deforms the neighbourhood $B_{r}(0)$ (a ball or radius $r$) of the origin on $T_{u}\Omega$ into a neighbourhood of $u$. Conversely, one can also regard the exponential map as an intrinsic local deformation (‘flattening’) of the manifold into the tangent space.

Geodesic distances

A result known as the Hopf-Rinow Theorem guarantees that geodesically complete manifolds are also complete metric spaces, in which one can realise a distance (called the geodesic distance or metric) between any pair of points $u,v$ as the length of the shortest path between them

$$d_{g}(u,v)=\min_{\gamma}\ell(\gamma)\quad\quad\text{s.t.}\quad\quad\gamma(0)=u,\,\,\gamma(T)=v,$$

which exists (i.e., the minimum is attained).

Isometries

Consider now a deformation of our manifold $\Omega$ into another manifold $\tilde{\Omega}$ with a Riemannian metric $h$, which we assume to be a diffeomorphism $\eta:(\Omega,g)\rightarrow(\tilde{\Omega},h)$ between the manifolds. Its differential $\mathrm{d}\eta:T\Omega\rightarrow T\tilde{\Omega}$ defines a map between the respective tangent bundles (referred to as pushforward), such that at a point $u$, we have $\mathrm{d}\eta_{u}:T_{u}\Omega\rightarrow T_{\eta(u)}\tilde{\Omega}$, interpreted as before: if we make a small displacement from point $u$ by tangent vector $X\in T_{u}\Omega$, the map $\eta$ will be displaced from point $\eta(u)$ by tangent vector $\mathrm{d}\eta_{u}(X)\in T_{\eta(u)}\tilde{\Omega}$.

Since the pushforwardprovides a mechanism to associate tangent vectors on the two manifolds, it allows to pullback the metric $h$ from $\tilde{\Omega}$ to $\Omega$,

$$(\eta^{*}h)_{u}(X,Y)=h_{\eta(u)}(\mathrm{d}\eta_{u}(X),\mathrm{d}\eta_{u}(Y))$$

If the pullback metric coincides at every point with that of $\Omega$, i.e., $g=\eta^{*}h$, the map $\eta$ is called (a Riemannian) isometry. For two-dimensional manifolds (surfaces), isometries can be intuitively understood as inelastic deformations that deform the manifold without ‘stretching’ or ‘tearing’ it.

By virtue of their definition, isometries preserve intrinsic structures such as geodesic distances, which are expressed entirely in terms of the Riemannian metric. Therefore, we can also understand isometries from the position of metric geometry, as distance-preserving maps (‘metric isometries’) between metric spaces $\eta:(\Omega,d_{g})\rightarrow(\tilde{\Omega},d_{h})$, in the sense that

$$d_{g}(u,v)=d_{h}(\eta(u),\eta(v))$$

for all $u,v\in\Omega$, or more compactly, $d_{g}=d_{h}\circ(\eta\times\eta)$. In other words, Riemannian isometries are also metric isometries. On connected manifolds, the converse is also true: every metric isometry is also a Riemannian isometry.

In our Geometric Deep Learning blueprint, $\eta$ is a model of domain deformations. When $\eta$ is an isometry, any intrinsic quantities are unaffected by such deformations. One can generalise exact (metric) isometries through the notions of metric dilation

$$\mathrm{dil}(\eta)=\sup_{u\neq v\in\Omega}\frac{d_{h}(\eta(u),\eta(v))}{d_{g}(u,v)}$$

or metric distortion

$$\mathrm{dis}(\eta)=\sup_{u,v\in\Omega}|d_{h}(\eta(u),\eta(v))-d_{g}(u,v)|,$$

which capture the relative and absolute change of the geodesic distances under $\eta$, respectively. The condition (5) for the stability of a function $f\in\mathcal{F}(\mathcal{X}(\Omega))$ under domain deformation can be rewritten in this case as

$$\|f(x,\Omega)-f(x\circ\eta^{-1},\tilde{\Omega})\|\leq C\|x\|\mathrm{dis}(\eta).$$

Intrinsic symmetries

A particular case of the above is a diffeomorphism of the domain itself (what we termed automorphism in Section 3.2), which we will denote by $\tau\in\mathrm{Diff}(\Omega)$. We will call it a Riemannian (self-)isometry if the pullback metric satisfies $\tau^{*}g=g$, or a metric (self-)isometry if $d_{g}=d_{g}\circ(\tau\times\tau)$. Not surprisingly,isometries form a group with the composition operator denoted by $\mathrm{Iso}(\Omega)$ and called the isometry group; the identity element is the map $\tau(u)=u$ and the inverse always exists (by definition of $\tau$ as a diffeomorphism). Self-isometries are thus intrinsic symmetries of manifolds.

Fourier analysis on Manifolds

We will now show how to construct intrinsic convolution-like operations on manifolds, which, by construction, will be invariant to isometric deformations. For this purpose, we have two options: One is to use an analogy of the Fourier transform, and define the convolution as a product in the Fourier domain. The other is to define the convolution spatially, by correlating a filter locally with the signal. Let us discuss the spectral approach first.

We remind that in the Euclidean domain the Fourier transform is obtained as the eigenvectors of circulant matrices, which are jointly diagonalisable due to their commutativity. Thus, any circulant matrix and in particular, differential operator, can be used to define an analogy of the Fourier transform on general domains. In Riemannian geometry, it is common to use the orthogonal eigenbasis of the Laplacian operator, which we will define here.

For this purpose, recall our definition of the intrinsic gradient operator $\nabla:\mathcal{X}(\Omega,\mathbb{R})\rightarrow\mathcal{X}(\Omega,T\Omega)$, producing a tangent vector field that indicates the local direction of steepest increase of a scalar field on the manifold. In a similar manner, we can define the divergence operator $\nabla^{*}:\mathcal{X}(\Omega,T\Omega)\rightarrow\mathcal{X}(\Omega,\mathbb{R})$. If we think of a tangent vector field as a flow on the manifold, the divergence measures the net flow of a field at a point, allowing to distinguish between field ‘sources’ and ‘sinks’. We use the notation $\nabla^{*}$ (as opposed to the common $\mathrm{div}$) to emphasise that the two operators are adjoint,

$$\langle X,\nabla x\rangle=\langle\nabla^{*}X,x\rangle,$$

where we use the inner products (15) and (16) between scalar and vector fields.

The Laplacian (also known as the Laplace-Beltrami operator in differential geometry) is an operator on $\mathcal{X}(\Omega)$ defined as $\Delta=\nabla^{*}\nabla$, which can be interpreted as the difference between the average of a function on an infinitesimal sphere around a point and the value of the function at the point itself. It is one of the most important operators in mathematical physics, used to describe phenomena as diverse as heat diffusion, quantum oscillations, and wave propagation. Importantly in our context, the Laplacian is intrinsic, and thus invariant under isometries of $\Omega$.

It is easy to see that the Laplacian is self-adjoint (‘symmetric’),

$$\langle\nabla x,\nabla x\rangle=\langle x,\Delta x\rangle=\langle\Delta x,x\rangle.$$

The quadratic form on the left in the above expression is actually the already familiar Dirichlet energy,

$$c^{2}(x)=\|\nabla x\|^{2}=\langle\nabla x,\nabla x\rangle=\int_{\Omega}\|\nabla x(u)\|_{u}^{2}\mathrm{d}u=\int_{\Omega}g_{u}(\nabla x(u),\nabla x(u))\mathrm{d}u$$

measuring the smoothness of $x$.

The Laplacian operator admits an eigedecomposition

$$\Delta\varphi_{k}=\lambda_{k}\varphi_{k},\quad\quad k=0,1,\ldots$$

with countable spectrum if the manifold is compact (which we tacitly assume), and orthogonal eigenfunctions, $\langle\varphi_{k},\varphi_{l}\rangle=\delta_{kl}$, due to the self-adjointness of $\Delta$. The Laplacian eigenbasis can also be constructed as a set of orthogonal minimisers of the Dirichlet energy,

$$\varphi_{k+1}=\arg\min_{\varphi}\|\nabla\varphi\|^{2}\quad\quad\text{s.t.}\quad\quad\|\varphi\|=1\,\,\,\text{and}\,\,\,\langle\varphi,\varphi_{j}\rangle=0$$

for $j=0,\ldots,k$, allowing to interpret it as the smoothest orthogonal basis on $\Omega$. The eigenfunctions $\varphi_{0},\varphi_{1},\ldots$ and the corresponding eigenvalues $0=\lambda_{0}\leq\lambda_{1}\leq\ldots$ can be interpreted as the analogy of the atoms and frequencies in the classical Fourier transform.

This orthogonal basis allows to expand square-integrable functions on $\Omega$ into Fourier series

$$x(u)=\sum_{k\geq 0}\langle x,\varphi_{k}\rangle\varphi_{k}(u)$$

where $\hat{x}_{k}=\langle x,\varphi_{k}\rangle$ are referred to as the Fourier coefficient or the (generalised) Fourier transform of $x$. Truncating the Fourier series results in an approximation error that can be bounded (Aflalo and Kimmel, 2013) by

$$\left\|x-\sum_{k=0}^{N}\langle x,\varphi_{k}\rangle\varphi_{k}\right\|^{2}\leq\frac{\|\nabla x\|^{2}}{\lambda_{N+1}}.$$

Aflalo et al. (2015) further showed that no other basis attains a better error, making the Laplacian eigenbasis optimal for representing smooth signals on manifolds.

Spectral Convolution on Manifolds

Spectral convolution can be defined as the product of Fourier transforms of the signal $x$ and the filter $\theta$,

$$(x\star\theta)(u)=\sum_{k\geq 0}(\hat{x}_{k}\cdot\hat{\theta}_{k})\varphi_{k}(u).$$ (17)

Note that here we use what is a property of the classical Fourier transform (the Convolution Theorem) as a way to define a non-Euclidean convolution. By virtue of its construction, the spectral convolution is intrinsic and thus isometry-invariant. Furthermore, since the Laplacian operator is isotropic, it has no sense of direction; in this sense, the situation is similar to that we had on graphs in Section 4.1 due to permutation invariance of neighbour aggregation.

Figure 12: Instability of spectral filters under domain perturbation. Left: a signal $\mathbf{x}$ on the mesh $\Omega$. Middle: result of spectral filtering in the eigenbasis of the Laplacian $\Delta$ on $\Omega$. Right: the same spectral filter applied to the eigenvectors of the Laplacian $\tilde{\Delta}$ of a nearly-isometrically perturbed domain $\tilde{\Omega}$ produces a very different result.

In practice, a direct computation of (17) appears to be prohibitively expensive due to the need to diagonalise the Laplacian. Even worse, it turns out geometrically unstable: the higher-frequency eigenfunctions of the Laplacian can change dramatically as a result of even small near-isometric perturbations of the domain $\Omega$ (see Figure 12). A more stable solution is provided by realising the filter as a spectral transfer function of the form $\hat{p}(\Delta)$,

$$(\hat{p}(\Delta)x)(u) = \sum_{k\geq 0}\hat{p}(\lambda_{k})\langle x,\varphi_{k}\rangle\varphi_{k}(u)$$ (18)

$$= \int_{\Omega}x(v)\,\sum_{k\geq 0}\hat{p}(\lambda_{k})\varphi_{k}(v)\varphi_{k}(u)\,\mathrm{d}v$$ (19)

which can be interpreted in two manners: either as a spectral filter (18), where we identify $\hat{\theta}_{k}=\hat{p}(\lambda_{k})$, or as a spatial filter (19) with a position-dependent kernel $\theta(u,v)=\sum_{k\geq 0}\hat{p}(\lambda_{k})\varphi_{k}(v)\varphi_{k}(u)$. The advantage of this formulation is that $\hat{p}(\lambda)$ can be parametrised by a small number of coefficients, and choosing parametric functions such as polynomials$\hat{p}(\lambda)=\sum_{l=0}^{r}\alpha_{l}\lambda^{l}$ allows for efficiently computing the filter as

$$(\hat{p}(\Delta)x)(u)=\sum_{k\geq 0}\sum_{l=0}^{r}\alpha_{l}\lambda_{k}^{l}\,\langle x,\varphi_{k}\rangle\varphi_{k}(u)=\sum_{l=0}^{r}\alpha_{l}(\Delta^{l}x)(u),$$

avoiding the spectral decomposition altogether. We will discuss this construction in further detail in Section 4.6.

Spatial Convolution on Manifolds

A second alternative is to attempt defining convolution on manifolds is by matching a filter at different points, like we did in formula (14),

$$(x\star\theta)(u)=\int_{T_{u}\Omega}x(\exp_{u}Y)\theta_{u}(Y)\mathrm{d}Y,$$ (20)

where we now have to use the exponential map to access the values of the scalar field $x$ from the tangent space, and the filter $\theta_{u}$ is defined in the tangent space at each point and hence position-dependent. If one defines the filter intrinsically, such a convolution would be isometry-invariant, a property we mentioned as crucial in many computer vision and graphics applications.

We need, however, to note several substantial differences from our previous construction in Sections 4.2–4.3. First, because a manifold is generally not a homogeneous space, we do not have anymore a global group structure allowing us have a shared filter (i.e., the same $\theta$ at every $u$ rather than $\theta_{u}$ in expression (20)) defined at one point and then move it around. An analogy of this operation on the manifold would require parallel transport, allowing to apply a shared $\theta$, defined as a function on $T_{u}\Omega$, at some other $T_{v}\Omega$. However, as we have seen, this in general will depend on the path between $u$ and $v$, so the way we move the filter around matters.Third, since we can use the exponential map only locally, the filter must be local, with support bounded by the injectivity radius. Fourth and most crucially, we cannot work with $\theta(X)$, as $X$ is an abstract geometric object: in order for it to be used for computations, we must represent it relative to some local basis $\omega_{u}:\mathbb{R}^{s}\rightarrow T_{u}\Omega$, as an $s$-dimensional array of coordinates $\mathbf{x}=\omega^{-1}_{u}(X)$. This allows us to rewrite the convolution (20) as

$$(x\star\theta)(u)=\int_{[0,1]^{s}}x(\exp_{u}(\omega_{u}\mathbf{y}))\theta(\mathbf{y})\mathrm{d}\mathbf{y},$$ (21)

with the filter defined on the unit cube. Since the exponential map is intrinsic (through the definition of geodesic), the resulting convolution is isometry-invariant.

Yet, this tacitly assumed we can carry the frame $\omega_{u}$ along to another manifold, i.e. $\omega_{u}^{\prime}=\mathrm{d}\eta_{u}\circ\omega_{u}$. Obtaining such a frame (or gauge, in physics terminology) given only a the manifold $\Omega$ in a consistent manner is however fraught with difficulty. First, a smooth global gauge may not exist: this is the situation on manifolds that are not parallelisable,in which case one cannot define a smooth non-vanishing tangent vector field. Second, we do not have a canonical gauge on manifolds, so this choice is arbitrary; since our convolution depends on $\omega$, if one chose a different one, we would obtain different results.

We should note that this is a case where practice diverges from theory: in practice, it is possible to build frames that are mostly smooth, with a limited number of singularities, e.g. by taking the intrinsic gradient of some intrinsic scalar field on the manifold. Moreover, such constructions are stable, i.e., the frames constructed this way will be identical on isometric manifolds and similar on approximately isometric ones. Such approaches were in fact employed in the early works on deep learning on manifolds (Masci et al., 2015; Monti et al., 2017).

Nevertheless, this solution is not entirely satisfactory because near singularities, the filter orientation (being defined in a fixed manner relative to the gauge) will vary wildly, leading to a non-smooth feature map even if the input signal and filter are smooth. Moreover, there is no clear reason why a given direction at some point $u$ should be considered equivalent to another direction at an altogether different point $v$. Thus, despite practical alternatives, we will look next for a more theoretically well-founded approach that would be altogether independent on the choice of gauge.

4.5 Gauges and Bundles

The notion of gauge, which we have defined as a frame for the tangent space, is quite a bit more general in physics: it can refer to a frame for anyvector bundle, not just the tangent bundle. Informally, a vector bundle describes a family of vector spaces parametrised by another space and consists of a base space $\Omega$ with an identical vector space $\mathbb{V}$ (called the fibre) attached to each position $u\in\Omega$ (for the tangent bundle these are the tangent spaces $T_{u}\Omega$). Roughly speaking, a bundle looks as a product $\Omega\times\mathbb{V}$ locally around $u$, but globally might be ‘twisted’ and have an overall different structure. In Geometric Deep Learning, fibres can be used to model the feature spaces at each point in the manifold $\Omega$, with the dimension of the fibre being equal to the number of feature channels. In this context, a new and fascinating kind of symmetry, called gauge symmetry may present itself.

Let us consider again an $s$-dimensional manifold $\Omega$ with its tangent bundle $T\Omega$, and a vector field $X:\Omega\rightarrow T\Omega$ (which in this terminology is referred to as a section on the tangent bundle). Relative to a gauge $\omega$ for the tangent bundle, $X$ is represented as a function $\mathbf{x}:\Omega\rightarrow\mathbb{R}^{s}$. However it is important to realise that what we are really interested in is the underlying geometrical object (vector field), whose representation as a function $\mathbf{x}\in\mathcal{X}(\Omega,\mathbb{R}^{s})$ depends on the choice of gauge $\omega$. If we change the gauge, we also need to change $\mathbf{x}$ so as to preserve the underlying vector field being represented.

Tangent bundles and the Structure group

When we change the gauge, we need to apply at each point an invertible matrix that maps the old gauge to the new one. This matrix is unique for every pair of gauges at each point, but possibly different at different points. In other words, a gauge transformation is a mapping ${\mathfrak{g}}:\Omega\rightarrow\operatorname{GL}(s)$, where $\operatorname{GL}(s)$ is the general linear group of invertible $s\times s$ matrices. It acts on the gauge $\omega_{u}:\mathbb{R}^{s}\rightarrow T_{u}\Omega$ to produce a new gauge $\omega^{\prime}_{u}=\omega_{u}\circ{\mathfrak{g}}_{u}:\mathbb{R}^{s}\rightarrow T_{u}\Omega$. The gauge transformation acts on a coordinate vector field at each point via $\mathbf{x}^{\prime}(u)={\mathfrak{g}}^{-1}_{u}\mathbf{x}(u)$ to produce the coordinate representation $\mathbf{x}^{\prime}$ of $X$ relative to the new gauge. The underlying vector field remains unchanged:

$$X(u)=\omega^{\prime}_{u}(\mathbf{x^{\prime}}(u))=\omega_{u}({\mathfrak{g}}_{u}{\mathfrak{g}}^{-1}_{u}\mathbf{x}(u))=\omega_{u}(\mathbf{x}(u))=X(u),$$

which is exactly the property we desired. More generally, we may have a field of geometric quantities that transform according to a representation $\rho$ of $\operatorname{GL}(s)$, e.g. a field of 2-tensors (matrices) $\mathbf{A}(u)\in\mathbb{R}^{s\times s}$ that transform like $\mathbf{A^{\prime}}(u)=\rho_{2}({\mathfrak{g}}^{-1}_{u})\mathbf{A}(u)=\rho_{1}({\mathfrak{g}}_{u})\mathbf{A}(u)\rho_{1}({\mathfrak{g}}^{-1}_{u})$. In this case, the gauge transformation ${\mathfrak{g}}_{u}$ acts via $\rho({\mathfrak{g}}_{u})$.

Sometimes we may wish to restrict attention to frames with a certain property, such as orthogonal frames, right-handed frames, etc. Unsurprisingly, we are interested in a set of some property-preserving transformations that form a group. For instance, the group that preserves orthogonality is the orthogonal group $\operatorname{O}(s)$ (rotations and reflections), and the group that additionally preserves orientation or ‘handedness’ is $\operatorname{SO}(s)$ (pure rotations). Thus, in general we have a group ${\mathfrak{G}}$ called the structure group of the bundle, and a gauge transformation is a map ${\mathfrak{g}}:\Omega\rightarrow{\mathfrak{G}}$. A key observation is that in all cases with the given property, for any two frames at a given point there exists exactly one gauge transformation relating them.

As mentioned before, gauge theory extends beyond tangent bundles, and in general, we can consider a bundle of vector spaces whose structure and dimensions are not necessarily related to those of the base space $\Omega$. For instance, a color image pixel has a position $u\in\Omega=\mathbb{Z}^{2}$ on a 2D grid and a value $\mathbf{x}(u)\in\mathbb{R}^{3}$ in the RGB space, so the space of pixels can be viewed as a vector bundle with base space $\mathbb{Z}^{2}$ and a fibre $\mathbb{R}^{3}$ attached at each point. It is customary to express an RGB image relative to a gauge that has basis vectors for R, G, and B (in that order), so that the coordinate representation of the image looks like $\mathbf{x}(u)=(r(u),g(u),b(u))^{\top}$. But we may equally well permute the basis vectors (color channels) independently at each position, as long as we remember the frame (order of channels) in use at each point. As a computational operation this is rather pointless, but as we will see shortly it is conceptually useful to think about gauge transformations for the space of RGB colors, because it allows us to express a gauge symmetry – in this case, an equivalence between colors – and make functions defined on images respect this symmetry (treating each color equivalently).

As in the case of a vector field on a manifold, an RGB gauge transformation changes the numerical representation of an image (permuting the RGB values independently at each pixel) but not the underlying image. In machine learning applications, we are interested in constructing functions $f\in\mathcal{F}(\mathcal{X}(\Omega))$ on such images (e.g. to perform image classification or segmentation), implemented as layers of a neural network. It follows that if, for whatever reason, we were to apply a gauge transformation to our image, we would need to also change the function $f$ (network layers) so as to preserve their meaning. Consider for simplicity a $1\times 1$ convolution, i.e. a map that takes an RGB pixel $\mathbf{x}(u)\in\mathbb{R}^{3}$ to a feature vector $\mathbf{y}(u)\in\mathbb{R}^{C}$. According to our Geometric Deep Learning blueprint, the output is associated with a group representation $\rho_{\textup{out}}$, in this case a $C$-dimensional representation of the structure group ${\mathfrak{G}}=\Sigma_{3}$ (RGB channel permutations), and similarly the input is associated with a representation $\rho_{\textup{in}}({\mathfrak{g}})={\mathfrak{g}}$. Then, if we apply a gauge transformation to the input, we would need to change the linear map ($1\times 1$ convolution) $f:\mathbb{R}^{3}\rightarrow\mathbb{R}^{C}$ to $f^{\prime}=\rho^{-1}_{\textup{out}}({\mathfrak{g}})\circ f\circ\rho_{\textup{in}}({\mathfrak{g}})$ so that the output feature vector $\mathbf{y}(u)=f(\mathbf{x}(u))$ transforms like $\mathbf{y}^{\prime}(u)=\rho_{\textup{out}}({\mathfrak{g}}_{u})\mathbf{y}(u)$ at every point. Indeed we verify:

$$\mathbf{y}^{\prime}=f^{\prime}(\mathbf{x}^{\prime})=\rho^{-1}_{\textup{out}}({\mathfrak{g}})f(\rho_{\textup{in}}({\mathfrak{g}})\rho^{-1}_{\textup{in}}({\mathfrak{g}})\mathbf{x})=\rho^{-1}_{\textup{out}}({\mathfrak{g}})f(\mathbf{x}).$$

Gauge Symmetries

To say that we consider gauge transformations to be symmetries is to say that any two gauges related by a gauge transformation are to be considered equivalent. For instance, if we take ${\mathfrak{G}}=\operatorname{SO}(d)$, any two right-handed orthogonal frames are considered equivalent, because we can map any such frame to any other such frame by a rotation. In other words, there are no distinguished local directions such as “up” or “right”. Similarly, if ${\mathfrak{G}}=\operatorname{O}(d)$ (the orthogonal group), then any left and right handed orthogonal frame are considered equivalent. In this case, there is no preferred orientation either. In general, we can consider a group ${\mathfrak{G}}$ and a collection of frames at every point $u$ such that for any two of them there is a unique ${\mathfrak{g}}(u)\in{\mathfrak{G}}$ that maps one frame onto the other.

Regarding gauge transformations as symmetries in our Geometric Deep Learning blueprint, we are interested in making the functions $f$ acting on signals defined on $\Omega$ and expressed with respect to the gauge should equivariant to such transformation. Concretely, this means that if we apply a gauge transformation to the input, the output should undergo the same transformation (perhaps acting via a different representation of ${\mathfrak{G}}$). We noted before that when we change the gauge, the function $f$ should be changed as well, but for a gauge equivariant map this is not the case: changing the gauge leaves the mapping invariant. To see this, consider again the RGB color space example. The map $f:\mathbb{R}^{3}\rightarrow\mathbb{R}^{C}$ is equivariant if $f\circ\rho_{\textup{in}}({\mathfrak{g}})=\rho_{\textup{out}}({\mathfrak{g}})\circ f$, but in this case the gauge transformation applied to $f$ has no effect: $\rho_{\textup{out}}^{-1}({\mathfrak{g}})\circ f\circ\rho_{\textup{in}}({\mathfrak{g}})=f$. In other words, the coordinate expression of a gauge equivariant map is independent of the gauge, in the same way that in the case of graph, we applied the same function regardless of how the input nodes were permuted. However, unlike the case of graphs and other examples covered so far, gauge transformations act not on $\Omega$ but separately on each of the feature vectors $\mathbf{x}(u)$ by a transformation ${\mathfrak{g}}(u)\in{\mathfrak{G}}$ for each $u\in\Omega$.

Further considerations enter the picture when we look at filters on manifolds with a larger spatial support. Let us first consider an easy example of a mapping $f:\mathcal{X}(\Omega,\mathbb{R})\rightarrow\mathcal{X}(\Omega,\mathbb{R})$ from scalar fields to scalar fields on an $s$-dimensional manifold $\Omega$. Unlike vectors and other geometric quantities, scalars do not have an orientation, so a scalar field $x\in\mathcal{X}(\Omega,\mathbb{R})$ is invariant to gauge transformations (it transforms according to the trivial representation $\rho({\mathfrak{g}})=1$). Hence, any linear map from scalar fields to scalar fields is gauge equivariant (or invariant, which is the same in this case). For example, we could write $f$ similarly to (19), as a convolution-like operation with a position-dependent filter $\theta:\Omega\times\Omega\rightarrow\mathbb{R}$,

$$(x\star\theta)(u)=\int_{\Omega}\theta(u,v)x(v)\mathrm{d}v.$$ (22)

This implies that we have a potentially different filter $\theta_{u}=\theta(u,\cdot)$ at each point, i.e., no spatial weight sharing — which gauge symmetry alone does not provide.

Consider now a more interesting case of a mapping $f:\mathcal{X}(\Omega,T\Omega)\rightarrow\mathcal{X}(\Omega,T\Omega)$ from vector fields to vector fields. Relative to a gauge, the input and output vector fields $X,Y\in\mathcal{X}(\Omega,T\Omega)$ are vector-valued functions $\mathbf{x},\mathbf{y}\in\mathcal{X}(\Omega,\mathbb{R}^{s})$. A general linear map between such functions can be written using the same equation we used for scalars (22), only replacing the scalar kernel by a matrix-valued one $\bm{\Theta}:\Omega\times\Omega\rightarrow\mathbb{R}^{s\times s}$. The matrix $\bm{\Theta}(u,v)$ should map tangent vectors in $T_{v}\Omega$ to tangent vectors in $T_{u}\Omega$, but these points have different gauges that we may change arbitrarily and independently. That is, the filter would have to satisfy $\bm{\Theta}(u,v)=\rho^{-1}({\mathfrak{g}}(u))\bm{\Theta}(u,v)\rho({\mathfrak{g}}(v))$ for all $u,v\in\Omega$, where $\rho$ denotes the action of ${\mathfrak{G}}$ on vectors, given by an $s\times s$ rotation matrix. Since ${\mathfrak{g}}(u)$ and ${\mathfrak{g}}(v)$ can be chosen freely, this is an overly strong constraint on the filter.

A better approach is to first transport the vectors to a common tangent space by means of the connection, and then impose gauge equivariance w.r.t. a single gauge transformation at one point only. Instead of (22), we can then define the following map between vector fields,

$$(\mathbf{x}\star\bm{\Theta})(u)=\int_{\Omega}\bm{\Theta}(u,v)\rho({\mathfrak{g}}_{v\rightarrow u})\mathbf{x}(v)\mathrm{d}v,$$ (23)

where ${\mathfrak{g}}_{v\rightarrow u}\in{\mathfrak{G}}$ denotes the parallel transport from $v$ to $u$ along the geodesic connecting these two points; its representation $\rho({\mathfrak{g}}_{v\rightarrow u})$ is an $s\times s$ rotation matrix rotating the vector as it moves between the points. Note that this geodesic is assumed to be unique, which is true only locally and thus the filter must have a local support. Under a gauge transformation ${\mathfrak{g}}_{u}$, this element transforms as ${\mathfrak{g}}_{u\rightarrow v}\mapsto{\mathfrak{g}}_{u}^{-1}{\mathfrak{g}}_{u\rightarrow v}{\mathfrak{g}}_{v}$, and the field itself transforms as $\mathbf{x}(v)\mapsto\rho({\mathfrak{g}}_{v})\mathbf{x}(v)$. If the filter commutes with the structure group representation $\bm{\Theta}(u,v)\rho({\mathfrak{g}}_{u})=\rho({\mathfrak{g}}_{u})\bm{\Theta}(u,v)$, equation (23) defines a gauge-equivariant convolution, which transforms as

$$(\mathbf{x}^{\prime}\star\bm{\Theta})(u)=\rho^{-1}({\mathfrak{g}}_{u})(\mathbf{x}\star\bm{\Theta})(u).$$

under the aforementioned transformation.

4.6 Geometric graphs and Meshes

We will conclude our discussion of different geometric domains with geometric graphs (i.e., graphs that can be realised in some geometric space) and meshes. In our ‘5G’ of geometric domains, meshes fall somewhere between graphs and manifolds: in many regards, they are similar to graphs, but their additional structure allows to also treat them similarly to continuous objects. For this reason, we do not consider meshes as a standalone object in our scheme, and in fact, will emphasise that many of the constructions we derive in this section for meshes are directly applicable to general graphs as well.

As we already mentioned in Section 4.4, two-dimensional manifolds (surfaces) are a common way of modelling 3D objects (or, better said, the boundary surfaces of such objects). In computer graphics and vision applications, such surfaces are often discretised as triangular meshes, which can be roughly thought of as a piece-wise planar approximation of a surface obtained by gluing triangles together along their edges. Meshes are thus (undirected) graphs with additional structure: in addition to nodes and edges, a mesh $\mathcal{T}=(\mathcal{V},\mathcal{E},\mathcal{F})$ also have ordered triplets of nodes forming triangular faces $\mathcal{F}=\{(u,v,q):u,v,q\in\mathcal{V}\,\,\,\text{and}\,\,\,(u,v),(u,q),(q,v)\in\mathcal{E}\}$; the order of the nodes defines the face orientation.

It is further assumed that that each edge is shared by exactly two triangles, and the boundary of all triangles incident on each node forms a single loop of edges. This condition guarantees that 1-hop neighbourhoods around each node are disk-like and the mesh thus constitutes a discrete manifold – such meshes are referred to as manifold meshes. Similarly to Riemannian manifolds, we can define a metric on the mesh. In the simplest instance, it can be induced from the embedding of the mesh nodes $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ and expressed through the Euclidean length of the edges, $\ell_{uv}=\|\mathbf{x}_{u}-\mathbf{x}_{v}\|$. A metric defined in this way automatically satisfies properties such as the triangle inequality, i.e., expressions of the form $\ell_{uv}\leq\ell_{uq}+\ell_{vq}$ for any $(u,v,q)\in\mathcal{F}$ and any combination of edges. Any property that can be expressed solely in terms of $\ell$ is intrinsic, and any deformation of the mesh preserving $\ell$ is an isometry – these notions are already familiar to the reader from our discussion in Section 4.4.

Laplacian matrices

By analogy to our treatment of graphs, let us assume a (manifold) mesh with $n$ nodes, each associated with a $d$-dimensional feature vector, which we can arrange (assuming some arbitrary ordering) into an $n\times d$ matrix $\mathbf{X}$. The features can represent the geometric coordinates of the nodes as well as additional properties such as colors, normals, etc, or in specific applications such as chemistry where geometric graphs model molecules, properties such as the atomic number.

Let us first look at the spectral convolution (17) on meshes, which we remind the readers, arises from the Laplacian operator. Considering the mesh as a discretisation of an underlying continuous surface, we can discretise the Laplacian as

$$(\bm{\Delta}\mathbf{X})_{u}=\sum_{v\in\mathcal{N}_{u}}w_{uv}(\mathbf{x}_{u}-\mathbf{x}_{v}),$$ (24)

or in matrix-vector notation, as an $n\times n$ symmetric matrix $\bm{\Delta}=\mathbf{D}-\mathbf{W}$, where $\mathbf{D}=\mathrm{diag}(d_{1},\ldots,d_{n})$ is called the degree matrix and $d_{u}=\sum_{v}w_{uv}$ the degree of node $u$. It is easy to see that equation (24) performs local permutation-invariant aggregation of neighbour features $\phi(\mathbf{x}_{u},\mathbf{X}_{\mathcal{N}_{u}})=d_{u}\mathbf{x}_{u}-\sum_{v\in\mathcal{N}_{u}}w_{uv}\mathbf{x}_{v}$, and $\mathbf{F}(\mathbf{X})=\bm{\Delta}\mathbf{X}$ is in fact an instance of our general blueprint (13) for constructing permutation-equivariant functions on graphs.

Note that insofar there is nothing specific to meshes in our definition of Laplacian in (24); in fact, this construction is valid for arbitrary graphs as well, with edge weights identified with the adjacency matrix, $\mathbf{W}=\mathbf{A}$, i.e., $w_{uv}=1$if $(u,v)\in\mathcal{E}$ and zero otherwise. Laplacians constructed in this way are often called combinatorial, to reflect the fact that they merely capture the connectivity structure of the graph. For geometric graphs (which do not necessarily have the additional structure of meshes, but whose nodes do have spatial coordinates that induces a metric in the form of edge lengths), it is common to use weights inversely related to the metric, e.g. $w_{uv}\propto e^{-\ell_{uv}}$.

On meshes, we can exploit the additional structure afforded by the faces, and define the edge weights in equation (24) using the cotangent formula (Pinkall and Polthier, 1993; Meyer et al., 2003)

$$w_{uv}=\frac{\cot\angle_{uqv}+\cot\angle_{upv}}{2a_{u}}$$ (25)

where $\angle_{uqv}$ and $\angle_{upv}$ are the two angles in the triangles $(u,q,v)$ and $(u,p,v)$ opposite the shared edge $(u,v)$, and $a_{u}$ is the local area element, typically computed as the area of the polygon constructed upon the barycenters of the triangles $(u,p,q)$ sharing the node $u$ and given by $a_{u}=\frac{1}{3}\sum_{v,q:(u,v,q)\in\mathcal{F}}a_{uvq}$.

The cotangent Laplacian can be shown to have multiple convenient properties (see e.g. Wardetzky et al. (2007)): it is a positive-semidefinite matrix, $\bm{\Delta}\succcurlyeq 0$ and thus has non-negative eigenvalues $\lambda_{1}\leq\ldots\leq\lambda_{n}$ that can be regarded as an analogy of frequency, it is symmetric and thus has orthogonal eigenvectors, and it is local (i.e., the value of $(\bm{\Delta}\mathbf{X})_{u}$ depends only on 1-hop neighbours, $\mathcal{N}_{u}$). Perhaps the most important property is the convergence of the cotangent mesh Laplacian matrix $\bm{\Delta}$ to the continuous operator $\Delta$ when the mesh is infinitely refined (Wardetzky, 2008). Equation (25) constitutes thus an appropriate discretisationof the Laplacian operator defined on Riemannian manifolds in Section 4.4.

While one expects the Laplacian to be intrinsic, this is not very obvious from equation (25), and it takes some effort to express the cotangent weights entirely in terms of the discrete metric $\ell$ as

$$w_{uv}=\frac{-\ell^{2}_{uv}+\ell^{2}_{vq}+\ell^{2}_{uq}}{8a_{uvq}}+\frac{-\ell^{2}_{uv}+\ell^{2}_{vp}+\ell^{2}_{up}}{8a_{uvp}}$$

where the area of the triangles $a_{ijk}$ is given as

$$a_{uvq}=\sqrt{s_{uvq}(s_{uvq}-\ell_{uv})(s_{uvq}-\ell_{vq})(s_{uvq}-\ell_{uq})}$$

using Heron’s semiperimeter formula with $s_{uvq}=\frac{1}{2}(\ell_{uv}+\ell_{uq}+\ell_{vq})$. This endows the Laplacian (and any quantities associated with it, such as its eigenvectors and eigenvalues) with isometry invariance, a property for which it is so loved in geometry processing and computer graphics (see an excellent review by Wang and Solomon (2019)): any deformation of the mesh that does not affect the metric $\ell$ (does not ‘stretch’ or ‘squeeze’ the edges of the mesh) does not change the Laplacian.

Finally, as we already noticed,the definition of the Laplacian (25) is invariant to the permutation of nodes in $\mathcal{N}_{u}$, as it involves aggregation in the form of summation. While on general graphs this is a necessary evil due to the lack of canonical ordering of neighbours, on meshes we can order the 1-hop neighbours according to some orientation (e.g., clock-wise), and the only ambiguity is the selection of the first node. Thus, instead of any possible permutation we need to account for cyclic shifts (rotations), which intuitively corresponds to the ambiguity arising from $\mathrm{SO}(2)$ gauge transformations discussed in Section 4.5. For a fixed gauge, it is possible to define an anisotropic Laplacian that is sensitive to local directions and amounts to changing the metric or the weights $w_{uv}$. Constructions of this kind were used to design shape descriptors by Andreux et al. (2014); Boscaini et al. (2016b) and in early Geometric Deep Learning architectures on meshes by Boscaini et al. (2016a).

Spectral analysis on meshes

The orthogonal eigenvectors $\bm{\Phi}=(\bm{\varphi}_{1},\ldots,\bm{\varphi}_{n})$ diagonalising the Laplacian matrix ($\bm{\Delta}=\bm{\Phi}\bm{\Lambda}\bm{\Phi}^{\top}$, where $\bm{\Lambda}=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{n})$ is the diagonal matrix of Laplacian eigenvalues), are used as the non-Euclidean analogy of the Fourier basis, allowing to perform spectral convolution on the mesh as the product of the respective Fourier transforms,

$$\mathbf{X}\star\bm{\theta}=\bm{\Phi}\,\mathrm{diag}(\bm{\Phi}^{\top}\bm{\theta})(\bm{\Phi}^{\top}\mathbf{X})=\bm{\Phi}\,\mathrm{diag}(\hat{\bm{\theta}})\hat{\mathbf{X}},$$

where the filter $\hat{\bm{\theta}}$ is designed directly in the Fourier domain. Again, nothing in this formula is specific to meshes, and one can use the Laplacian matrix of a generic (undirected) graph.It is tempting to exploit this spectral definition of convolution to generalise CNNs to graphs, which in fact was done by one of the authors of this text, Bruna et al. (2013). However, it appears that the non-Euclidean Fourier transform is extremely sensitive to even minor perturbations of the underlying mesh or graph (see Figure 12 in Section 4.4) and thus can only be used when one has to deal with different signals on a fixed domain, but not when one wishes to generalise across different domains. Unluckily, many computer graphics and vision problems fall into the latter category, where one trains a neural network on one set of 3D shapes (meshes) and test on a different set, making the Fourier transform-based approach inappropriate.

As noted in Section 4.4, it is preferable to use spectral filters of the form (18) applying some transfer function $\hat{p}(\lambda)$ to the Laplacian matrix,

$$\hat{p}(\bm{\Delta})\mathbf{X}=\bm{\Phi}\hat{p}(\bm{\Lambda})\bm{\Phi}^{\top}\mathbf{X}=\bm{\Phi}\,\mathrm{diag}(\hat{p}(\lambda_{1}),\ldots,\hat{p}(\lambda_{n}))\hat{\mathbf{X}}.$$

When $\hat{p}$ can be expressed in terms of matrix-vector products, the eigendecomposition of the $n\times n$ matrix $\bm{\Delta}$ can be avoided altogether. For example, Defferrard et al. (2016) used polynomials of degree $r$ as filter functions,

$$\hat{p}(\bm{\Delta})\mathbf{X}=\sum_{k=0}^{r}\alpha_{k}\bm{\Delta}^{k}\mathbf{X}=\alpha_{0}\mathbf{X}+\alpha_{1}\bm{\Delta}\mathbf{X}+\ldots+\alpha_{r}\bm{\Delta}^{r}\mathbf{X},$$

amounting to the multiplication of the $n\times d$ feature matrix $\mathbf{X}$ by the $n\times n$ Laplacian matrix $r$ times. Since the Laplacian is typically sparse (with $\mathcal{O}(|\mathcal{E}|)$ non-zero elements) this operation has low complexity of $\mathcal{O}(|\mathcal{E}|dr)\sim\mathcal{O}(|\mathcal{E}|)$. Furthermore, since the Laplacian is local, a polynomial filter of degree $r$ is localised in $r$-hop neighbourhood.

However, this exact property comes at a disadvantage when dealing with meshes, since the actual support of the filter (i.e., the radius it covers) depends on the resolution of the mesh. One has to bear in mind that meshes arise from the discretisation of some underlying continuous surface, and one may have two different meshes $\mathcal{T}$ and $\mathcal{T}^{\prime}$ representing the same object. In a finer mesh, one might have to use larger neighbourhoods (thus, larger degree $r$ of the filter) than in a coarser one.

For this reason, in computer graphics applications it is more common to use rational filters, since they are resolution-independent. There are many ways to define such filters (see, e.g. Patanè (2020)), the most common being as a polynomial of some rational function, e.g., $\frac{\lambda-1}{\lambda+1}$. More generally, one can use a complex function, such as the Cayley transform $\frac{\lambda-\mathrm{i}}{\lambda+\mathrm{i}}$ that maps the real line into the unit circle in the complex plane. Levie et al. (2018) used spectral filters expressed as Cayley polynomials, real rational functions with complex coefficients $\alpha_{l}\in\mathbb{C}$,

$$\hat{p}(\lambda)=\mathrm{Re}\left(\sum_{l=0}^{r}\alpha_{l}\left(\frac{\lambda-\mathrm{i}}{\lambda+\mathrm{i}}\right)^{l}\right).$$

When applied to matrices, the computation of the Cayley polynomial requires matrix inversion,

$$\hat{p}(\bm{\Delta})=\mathrm{Re}\left(\sum_{l=0}^{r}\alpha_{l}(\bm{\Delta}-\mathrm{i}\mathbf{I})^{l}(\bm{\Delta}+\mathrm{i}\mathbf{I})^{-l}\right),$$

which can be carried out approximately with linear complexity. Unlike polynomial filters, rational filters do not have a local support, but have exponential decay (Levie et al., 2018). A crucial difference compared to the direct computation of the Fourier transform is that polynomial and rational filters are stable under approximate isometric deformations of the underlying graph or mesh – various results of this kind were shown e.g. by Levie et al. (2018, 2019); Gama et al. (2020); Kenlay et al. (2021).

Meshes as operators and Functional maps

The paradigm of functional maps suggests thinking of meshes as operators. As we will show, this allows obtaining more interesting types of invariance exploiting the additional structure of meshes. For the purpose of our discussion, assume the mesh $\mathcal{T}$ is constructed upon embedded nodes with coordinates $\mathbf{X}$. If we construct an intrinsic operator like the Laplacian, it can be shown that it encodes completely the structure of the mesh, and one can recover the mesh (up to its isometric embedding, as shown by Zeng et al. (2012)). This is also true for some other operators (see e.g. Boscaini et al. (2015); Corman et al. (2017); Chern et al. (2018)), so we will assume a general operator, or $n\times n$ matrix $\mathbf{Q}(\mathcal{T},\mathbf{X})$, as a representation of our mesh.

In this view, the discussion of Section 4.1 of learning functions of the form $f(\mathbf{X},\mathcal{T})$ can be rephrased as learning functions of the form $f(\mathbf{Q})$. Similar to graphs and sets, the nodes of meshes also have no canonical ordering, i.e., functions on meshes must satisfy the permutation invariance or equivariance conditions,

$$f(\mathbf{Q}) = f(\mathbf{P}\mathbf{Q}\mathbf{P}^{\top})$$

$$\mathbf{P}\mathbf{F}(\mathbf{Q}) = \mathbf{F}(\mathbf{P}\mathbf{Q}\mathbf{P}^{\top})$$

for any permutation matrix $\mathbf{P}$. However, compared to general graphs we now have more structure: we can assume that our mesh arises from the discretisation of some underlying continuous surface $\Omega$. It is thus possible to have a different mesh $\mathcal{T}^{\prime}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime},\mathcal{F}^{\prime})$ with $n^{\prime}$ nodes and coordinates $\mathbf{X}^{\prime}$ representing the same object $\Omega$ as $\mathcal{T}$. Importantly, the meshes $\mathcal{T}$ and $\mathcal{T}^{\prime}$ can have a different connectivity structure and even different number of nodes ($n^{\prime}\neq n$). Therefore, we cannot think of these meshes as isomorphic graphs with mere reordering of nodes and consider the permutation matrix $\mathbf{P}$ as correspondence between them.

Functional maps were introduced by Ovsjanikov et al. (2012) as a generalisation of the notion of correspondence to such settings, replacing the correspondence between points on two domains (a map $\eta:\Omega\rightarrow\Omega^{\prime}$) with correspondence between functions (a map $\mathbf{C}:\mathcal{X}(\Omega)\rightarrow\mathcal{X}(\Omega^{\prime})$, see Figure 13). A functional map is a linear operator $\mathbf{C}$, represented as a matrix $n^{\prime}\times n$, establishing correspondence between signals $\mathbf{x}^{\prime}$ and $\mathbf{x}$ on the respective domains as

$$\mathbf{x}^{\prime}=\mathbf{C}\mathbf{x}.$$

Rustamov et al. (2013) showed that in order to guarantee area-preserving mapping, the functional map must be orthogonal, $\mathbf{C}^{\top}\mathbf{C}=\mathbf{I}$, i.e., be an element of the orthogonal group $\mathbf{C}\in\mathrm{O}(n)$. In this case, we can invert the map using $\mathbf{C}^{-1}=\mathbf{C}^{\top}$.

Figure 13: Pointwise map (left) vs functional map (right).

The functional map also establishes a relation between the operator representation of meshes,

$$\mathbf{Q}^{\prime}=\mathbf{C}\mathbf{Q}\mathbf{C}^{\top},\quad\quad\mathbf{Q}=\mathbf{C}^{\top}\mathbf{Q}^{\prime}\mathbf{C},$$

which we can interpret as follows: given an operator representation $\mathbf{Q}$ of $\mathcal{T}$ and a functional map $\mathbf{C}$, we can construct its representation $\mathbf{Q}^{\prime}$ of $\mathcal{T}^{\prime}$ by first mapping the signal from $\mathcal{T}^{\prime}$ to $\mathcal{T}$ (using $\mathbf{C}^{\top}$), applying the operator $\mathbf{Q}$, and then mapping back to $\mathcal{T}^{\prime}$ (using $\mathbf{C}$)This leads us to a more general class of remeshing invariant (or equivariant) functions on meshes, satisfying

$$f(\mathbf{Q}) = f(\mathbf{C}\mathbf{Q}\mathbf{C}^{\top})=f(\mathbf{Q}^{\prime})$$

$$\mathbf{C}\mathbf{F}(\mathbf{Q}) = \mathbf{F}(\mathbf{C}\mathbf{Q}\mathbf{C}^{\top})=\mathbf{F}(\mathbf{Q}^{\prime})$$

for any $\mathbf{C}\in\mathrm{O}(n)$. It is easy to see that the previous setting of permutation invariance and equivariance is a particular case, which can be thought of as a trivial remeshing in which only the order of nodes is changed.

Wang et al. (2019a) showed that given an eigendecomposition of the operator $\mathbf{Q}=\mathbf{V}\bm{\Lambda}\mathbf{V}^{\top}$, any remeshing invariant (or equivariant) function can be expressed as $f(\mathbf{Q})=f(\bm{\Lambda})$ and $\mathbf{F}(\mathbf{Q})=\mathbf{V}\mathbf{F}(\bm{\Lambda})$, or in other words, remeshing-invariant functions involve only the spectrum of $\mathbf{Q}$. Indeed, functions of Laplacian eigenvalues have been proven in practice to be robust to surface discretisation and perturbation, explaining the popularity of spectral constructions based on Laplacians in computer graphics, as well as in deep learning on graph (Defferrard et al., 2016; Levie et al., 2018). Since this result refers to a generic operator $\mathbf{Q}$, multiple choices are available besides the ubiquitous Laplacian – notable examples include the Dirac (Liu et al., 2017; Kostrikov et al., 2018) or Steklov (Wang et al., 2018) operators, as well as learnable parametric operators (Wang et al., 2019a).

5 Geometric Deep Learning Models

Having thoroughly studied various instantiations of our Geometric Deep Learning blueprint (for different choices of domain, symmetry group, and notions of locality), we are ready to discuss how enforcing these prescriptions can yield some of the most popular deep learning architectures.

Our exposition, once again, will not be in strict order of generality. We initially cover three architectures for which the implementation follows nearly-directly from our preceding discussion: convolutional neural networks (CNNs), group-equivariant CNNs, and graph neural networks (GNNs).

We will then take a closer look into variants of GNNs for cases where a graph structure is not known upfront (i.e. unordered sets), and through our discussion we will describe the popular Deep Sets and Transformer architectures as instances of GNNs.

Following our discussion on geometric graphs and meshes, we first describe equivariant message passing networks, which introduce explicit geometric symmetries into GNN computations. Then, we show ways in which our theory of geodesics and gauge symmetries can be materialised within deep learning, recovering a family of intrinsic mesh CNNs (including Geodesic CNNs, MoNet and gauge-equivariant mesh CNNs).

Lastly, we look back on the grid domain from a temporal angle. This discussion will lead us to recurrent neural networks (RNNs). We will demonstrate a manner in which RNNs are translation equivariant over temporal grids, but also study their stability to time warping transformations. This property is highly desirable for properly handling long-range dependencies, and enforcing class invariance to such transformations yields exactly the class of gated RNNs (including popular RNN models such as the LSTM or GRU).

While we hope the above canvasses most of the key deep learning architectures in use at the time of writing, we are well aware that novel neural network instances are proposed daily. Accordingly, rather than aiming to cover every possible architecture, we hope that the following sections are illustrative enough, to the point that the reader is able to easily categorise any future Geometric Deep Learning developments using the lens of invariances and symmetries.

5.1 Convolutional Neural Networks

Convolutional Neural Networks are perhaps the earliest and most well known example of deep learning architectures following the blueprint of Geometric Deep Learning outlined in Section 3.5. In Section 4.2 we have fully characterised the class of linear and local translation equivariant operators, given by convolutions ${\bf C(\bm{\theta})}{\bf x}={\bf x}\star\bm{\theta}$ with a localised filter $\bm{\theta}$. Let us first focus on scalar-valued (‘single-channel’ or ‘grayscale’) discretised images, where the domain is the grid $\Omega=[H]\times[W]$ with $\mathbf{u}=(u_{1},u_{2})$ and ${\bf x}\in{\mathcal{X}}(\Omega,\mathbb{R})$.

Any convolution with a compactly supported filter of size $H^{f}\times W^{f}$ can be written as a linear combination of generators $\bm{\theta}_{1,1},\dots,\bm{\theta}_{{H^{f},W^{f}}}$, given for example by the unit peaks $\bm{\theta}_{vw}(u_{1},u_{2})=\delta(u_{1}-v,u_{2}-w)$. Any local linear equivariant map is thus expressible as

$$\mathbf{F}(\mathbf{x})=\sum_{v=1}^{H^{f}}\sum_{w=1}^{W^{f}}\alpha_{vw}\mathbf{C}(\bm{\theta}_{vw})\mathbf{x}~{},$$ (26)

which, in coordinates, corresponds to the familiar 2D convolution (see Figure 14 for an overview):

Figure 14: The process of convolving an image $\mathbf{x}$ with a filter $\mathbf{C}(\bm{\theta})$. The filter parameters $\bm{\theta}$ can be expressed as a linear combination of generators $\bm{\theta}_{vw}$.

$$\mathbf{F}(\mathbf{x})_{uv}=\sum_{a=1}^{H^{f}}\sum_{b=1}^{W^{f}}{\alpha}_{ab}x_{u+a,v+b}~{}.$$ (27)

Other choices of the basis $\bm{\theta}_{vw}$ are also possible and will yield equivalent operations (for potentially different choices of $\alpha_{vw}$). A popular example are directional derivatives: $\bm{\theta}_{vw}(u_{1},u_{2})=\delta(u_{1},u_{2})-\delta(u_{1}-v,u_{2}-w),(v,w)\neq(0,0)$ taken together with the local average $\bm{\theta}_{0}(u_{1},u_{2})=\frac{1}{H_{f}W_{f}}$. In fact, directional derivatives can be considered a grid-specific analogue of diffusion processes on graphs, which we recover if we assume each pixel to be a node connected to its immediate neighbouring pixels in the grid.

When the scalar input channel is replaced by multiple channels (e.g., RGB colours, or more generally an arbitrary number of feature maps), the convolutional filter becomes a convolutional tensor expressing arbitrary linear combinations of input features into output feature maps. In coordinates, this can be expressed as

$$\mathbf{F}(\mathbf{x})_{uvj}=\sum_{a=1}^{H^{f}}\sum_{b=1}^{W^{f}}\sum_{c=1}^{M}{\alpha}_{jabc}x_{u+a,v+b,c}~{},~{}j\in[N]~{},$$ (28)

where $M$ and $N$ are respectively the number of input and output channels. This basic operation encompasses a broad class of neural network architectures, which, as we will show in the next section, have had a profound impact across many areas of computer vision, signal processing, and beyond. Here, rather than dissecting the myriad of possible architectural variants of CNNs, we prefer to focus on some of the essential innovations that enabled their widespread use.

Efficient multiscale computation

As discussed in the GDL template for general symmetries, extracting translation invariant features out of the convolutional operator $\mathbf{F}$ requires a non-linear step.Convolutional features are processed through a non-linear activation function $\sigma$, acting element-wise on the input—i.e., $\sigma:{\mathcal{X}}(\Omega)\to{\mathcal{X}}(\Omega)$, as $\sigma(\mathbf{x})(u)=\sigma(\mathbf{x}(u))$. Perhaps the most popular example at the time of writing is the Rectified Linear Unit (ReLU): $\sigma(x)=\max(x,0)$. This non-linearity effectively rectifies the signals, pushing their energy towards lower frequencies, and enabling the computation of high-order interactions across scales by iterating the construction.

Already in the early works of Fukushima and Miyake (1982) and LeCun et al. (1998), CNNs and similar architectures had a multiscale structure, where after each convolutional layer (28) one performs a grid coarsening $\mathbf{P}:{\mathcal{X}}(\Omega)\to{\mathcal{X}}(\Omega^{\prime})$, where the grid $\Omega^{\prime}$ has coarser resolution than $\Omega$. This enables multiscale filters with effectively increasing receptive field, yet retaining a constant number of parameters per scale. Several signal coarsening strategies $\mathbf{P}$ (referred to as pooling) may be used, the most common are applying a low-pass anti-aliasing filter (e.g. local average) followed by grid downsampling, or non-linear max-pooling.

In summary, a ‘vanilla’ CNN layer can be expressed as the composition of the basic objects already introduced in our Geometric Deep Learning blueprint:

$$\mathbf{h}=\mathbf{P}(\sigma(\mathbf{F}(\mathbf{x})))~{},$$ (29)

i.e. an equivariant linear layer $\mathbf{F}$, a coarsening operation $\mathbf{P}$, and a non-linearity $\sigma$. It is also possible to perform translation invariant global pooling operations within CNNs. Intuitively, this involves each pixel—which, after several convolutions, summarises a patch centered around it—proposing the final representation of the image, with the ultimate choice being guided by a form of aggregation of these proposals. A popular choice here is the average function, as its outputs will retain similar magnitudes irrespective of the image size (Springenberg et al., 2014).

Prominent examples following this CNN blueprint (some of which we will discuss next) are displayed in Figure 15.

Figure 15: Prominent examples of CNN architectures. Top-to-bottom: LeNet (LeCun et al., 1998), AlexNet (Krizhevsky et al., 2012), ResNet (He et al., 2016) and U-Net (Ronneberger et al., 2015). Drawn using the PlotNeuralNet package (Iqbal, 2018).

Deep and Residual Networks

A CNN architecture, in its simplest form, is therefore specified by hyperparameters $(H^{f}_{k},W^{f}_{k},N_{k},p_{k})_{k\leq K}$, with $M_{k+1}=N_{k}$ and $p_{k}=0,1$ indicating whether grid coarsening is performed or not. While all these hyperparameters are important in practice, a particularly important question is to understand the role of depth $K$ in CNN architectures, and what are the fundamental tradeoffs involved in choosing such a key hyperparameter, especially in relation to the filter sizes ($H_{k}^{f},W_{k}^{f}$).

While a rigorous answer to this question is still elusive, mounting empirical evidence collected throughout the recent years suggests a favourable tradeoff towards deeper (large $K$) yet thinner (small $(H^{f}_{k},W^{f}_{k})$) models . In this context, a crucial insight by He et al. (2016) was to reparametrise each convolutional layer to model a perturbation of the previous features, rather than a generic non-linear transformation:

$$\mathbf{h}=\mathbf{P}\left(\mathbf{x}+\sigma(\mathbf{F}(\mathbf{x}))\right)~{}.$$ (30)

The resulting residual networks provide several key advantages over the previous formulation. In essence, the residual parametrisation is consistent with the view that the deep network is a discretisation of an underlying continuous dynamical system, modelled as an ordinary differential equation (ODE). Crucially, learning a dynamical system by modeling its velocity turns out to be much easier than learning its position directly. In our learning setup, this translates into an optimisation landscape with more favorable geometry, leading to the ability to train much deeper architectures than was possible before. As will be discussed in future work, learning using deep neural networks defines a non-convex optimisation problem, which can be efficiently solved using gradient-descent methods under certain simplifying regimes. The key advantage of the ResNet parametrisation has been rigorously analysed in simple scenarios (Hardt and Ma, 2016), and remains an active area of theoretical investigation. Finally, Neural ODEs (Chen et al., 2018) are a recent popular architecture that pushes the analogy with ODEs even further, by learning the parameters of the ODE $\dot{\mathbf{x}}=\sigma(\mathbf{F}(\mathbf{x}))$ directly and relying on standard numerical integration.

Normalisation

Another important algorithmic innovation that boosted the empirical performance of CNNs significantly is the notion of normalisation. In early models of neural activity, it was hypothesised that neurons perform some form of local ‘gain control’, where the layer coefficients $\mathbf{x}_{k}$ are replaced by $\mathbf{\tilde{x}}_{k}=\sigma_{k}^{-1}\odot(\mathbf{x}_{k}-\mu_{k})$. Here, $\mu_{k}$ and $\sigma_{k}$ encode the first and second-order moment information of $\mathbf{x}_{k}$, respectively. Further, they can be either computed globally or locally.

In the context of Deep Learning, this principle was widely adopted through the batch normalisation layer (Ioffe and Szegedy, 2015), followed by several variants (Ba et al., 2016; Salimans and Kingma, 2016; Ulyanov et al., 2016; Cooijmans et al., 2016; Wu and He, 2018). Despite some attempts to rigorously explain the benefits of normalisation in terms of better conditioned optimisation landscapes (Santurkar et al., 2018), a general theory that can provide guiding principles is still missing at the time of writing.

Data augmentation

While CNNs encode the geometric priors associated with translation invariance and scale separation, they do not explicitly account for other known transformations that preserve semantic information, e.g lightning or color changes, or small rotations and dilations. A pragmatic approach to incorporate these priors with minimal architectural changes is to perform data augmentation, where one manually performs said transformations to the input images and adds them into the training set.

Data augmentation has been empirically successful and is widely used—not only to train state-of-the-art vision architectures, but also to prop up several developments in self-supervised and causal representation learning (Chen et al., 2020; Grill et al., 2020; Mitrovic et al., 2020). However, it is provably sub-optimal in terms of sample complexity (Mei et al., 2021); a more efficient strategy considers instead architectures with richer invariance groups—as we discuss next.

5.2 Group-equivariant CNNs

As discussed in Section 4.3, we can generalise the convolution operation from signals on a Euclidean space to signals on any homogeneous space $\Omega$ acted upon by a group ${\mathfrak{G}}$. By analogy to the Euclidean convolution where a translated filter is matched with the signal, the idea of group convolution is to move the filter around the domain using the group action, e.g. by rotating and translating. By virtue of the transitivity of the group action, we can move the filter to any position on $\Omega$. In this section, we will discuss several concrete examples of the general idea of group convolution, including implementation aspects and architectural choices.

Discrete group convolution

We begin by considering the case where the domain $\Omega$ as well as the group ${\mathfrak{G}}$ are discrete. As our first example, we consider medical volumetric images represented as signals of on 3D grids with discrete translation and rotation symmetries. The domain is the 3D cubical grid $\Omega=\mathbb{Z}^{3}$ and the images (e.g. MRI or CT 3D scans) are modelled as functions $x:\mathbb{Z}^{3}\rightarrow\mathbb{R}$, i.e. $x\in\mathcal{X}(\Omega)$. Although in practice such images have support on a finite cuboid $[W]\times[H]\times[D]\subset\mathbb{Z}^{3}$, we instead prefer to view them as functions on $\mathbb{Z}^{3}$ with appropriate zero padding. As our symmetry, we consider the group ${\mathfrak{G}}=\mathbb{Z}^{3}\rtimes O_{h}$ of distance- and orientation-preserving transformations on $\mathbb{Z}^{3}$. This group consists of translations ($\mathbb{Z}^{3}$) and the discrete rotations $O_{h}$ generated by $90$ degree rotations about the three axes (see Figure 16).

Figure 16: A $3\times 3$ filter, rotated by all $24$ elements of the discrete rotation group $O_{h}$, generated by $90$-degree rotations about the vertical axis (red arrows), and $120$-degree rotations about a diagonal axis (blue arrows).

As our second example, we consider DNAsequences made up of four letters: C, G, A, and T. The sequences can be represented on the 1D grid $\Omega=\mathbb{Z}$ as signals $x:\mathbb{Z}\rightarrow\mathbb{R}^{4}$, where each letter is one-hot coded in $\mathbb{R}^{4}$. Naturally, we have a discrete 1D translation symmetry on the grid, but DNA sequences have an additional interesting symmetry. This symmetry arises from the way DNA is physically embodied as a double helix, and the way it is read by the molecular machinery of the cell. Each strand of the double helix begins with what is called the $5^{\prime}$-end and ends with a $3^{\prime}$-end, with the $5^{\prime}$ on one strand complemented by a $3^{\prime}$ on the other strand. In other words, the two strands have an opposite orientation.Since the DNA molecule is always read off starting at the $5^{\prime}$-end, but we do not know which one, a sequence such as ACCCTGG is equivalent to the reversed sequence with each letter replaced by its complement, CCAGGGT. This is called reverse-complement symmetry of the letter sequence. We thus have the two-element group $\mathbb{Z}_{2}=\{0,1\}$ corresponding to the identity $0$ and reverse-complement transformation $1$ (and composition $1+1=0\mod{2}$). The full group combines translations and reverse-complement transformations.

In our case, the group convolution (14) we defined in Section 4.3 is given as

$$(x\star\theta)({\mathfrak{g}})=\sum_{u\in\Omega}x_{u}\rho({\mathfrak{g}})\theta_{u},$$ (31)

the inner product between the (single-channel) input signal $x$ and a filter $\theta$ transformed by ${\mathfrak{g}}\in{\mathfrak{G}}$ via $\rho({\mathfrak{g}})\theta_{u}=\theta_{{\mathfrak{g}}^{-1}u}$, and the output $x\star\theta$ is a function on ${\mathfrak{G}}$. Note that since $\Omega$ is discrete, we have replaced the integral from equation (14) by a sum.

Transform+Convolve approach

We will show that the group convolution can be implemented in two steps: a filter transformation step, and a translational convolution step. The filter transformation step consists of creating rotated (or reverse-complement transformed) copies of a basic filter, while the translational convolution is the same as in standard CNNs and thus efficiently computable on hardware such as GPUs. To see this, note that in both of our examples we can write a general transformation ${\mathfrak{g}}\in{\mathfrak{G}}$ as a transformation ${\mathfrak{h}}\in{\mathfrak{H}}$ (e.g. a rotation or reverse-complement transformation) followed by a translation ${\mathfrak{k}}\in\mathbb{Z}^{d}$, i.e. ${\mathfrak{g}}={\mathfrak{k}}{\mathfrak{h}}$ (with juxtaposition denoting the composition of the group elements ${\mathfrak{k}}$ and ${\mathfrak{h}}$). By properties of the group representation, we have $\rho({\mathfrak{g}})=\rho({\mathfrak{k}}{\mathfrak{h}})=\rho({\mathfrak{k}})\rho({\mathfrak{h}})$. Thus,

$$(x\star\theta)({\mathfrak{k}}{\mathfrak{h}}) =\sum_{u\in\Omega}x_{u}\rho({\mathfrak{k}})\rho({\mathfrak{h}})\theta_{u}$$ (32)

$$=\sum_{u\in\Omega}x_{u}(\rho({\mathfrak{h}})\theta)_{u-{\mathfrak{k}}}$$

We recognise the last equation as the standard (planar Euclidean) convolution of the signal $x$ and the transformed filter $\rho({\mathfrak{h}})\theta$. Thus, to implement group convolution for these groups, we take the canonical filter $\theta$, create transformed copies $\theta_{\mathfrak{h}}=\rho({\mathfrak{h}})\theta$ for each ${\mathfrak{h}}\in{\mathfrak{H}}$ (e.g. each rotation ${\mathfrak{h}}\in O_{h}$ or reverse-complement DNA symmetry ${\mathfrak{h}}\in\mathbb{Z}_{2}$), and then convolve $x$ with each of these filters: $(x\star\theta)({\mathfrak{k}}{\mathfrak{h}})=(x\star\theta_{{\mathfrak{h}}})({\mathfrak{k}})$. For both of our examples, the symmetries act on filters by simply permuting the filter coefficients, as shown in Figure 16 for discrete rotations. Hence, these operations can be implemented efficiently using an indexing operation with pre-computed indices.

While we defined the feature maps output by the group convolution $x\star\theta$ as functions on ${\mathfrak{G}}$, the fact that we can split ${\mathfrak{g}}$ into ${\mathfrak{h}}$ and ${\mathfrak{k}}$ means that we can also think of them as a stack of Euclidean feature maps (sometimes called orientation channels), with one feature map per filter transformation / orientation ${\mathfrak{k}}$. For instance, in our first example we would associate to each filter rotation (each node in Figure 16) a feature map, which is obtained by convolving (in the traditional translational sense) the rotated filter. These feature maps can thus still be stored as a $W\times H\times C$ array, where the number of channels $C$ equals the number of independent filters times the number of transformations ${\mathfrak{h}}\in{\mathfrak{H}}$ (e.g. rotations).

As shown in Section 4.3, the group convolution is equivariant: $(\rho({\mathfrak{g}})x)\star\theta=\rho({\mathfrak{g}})(x\star\theta)$. What this means in terms of orientation channels is that under the action of ${\mathfrak{h}}$, each orientation channel is transformed, and the orientation channels themselves are permuted. For instance, if we associate one orientation channel per transformation in Figure 16 and apply a rotation by $90$ degrees about the z-axis (corresponding to the red arrows), the feature maps will be permuted as shown by the red arrows. This description makes it clear that a group convolutional neural network bears much similarity to a traditional CNN. Hence, many of the network design patterns discussed in the Section 5.1, such as residual networks, can be used with group convolutions as well.

Spherical CNNs in the Fourier domain

For the continuous symmetry group of the sphere that we saw in Section 4.3, it is possible to implement the convolution in the spectral domain, using the appropriate Fourier transform (we remind the reader that the convolution on $\mathbb{S}^{2}$ is a function on $\mathrm{SO}(3)$, hence we need to define the Fourier transform on both these domains in order to implement multi-layer spherical CNNs). Spherical harmonics are an orthogonal basis on the 2D sphere, analogous to the classical Fourier basis of complex exponential. On the special orthogonal group, the Fourier basis is known as the Wigner D-functions. In both cases, the Fourier transforms (coefficients) are computed as the inner product with the basis functions, and an analogy of the Convolution Theorem holds: one can compute the convolution in the Fourier domain as the element-wise product of the Fourier transforms. Furthermore, FFT-like algorithms exist for the efficient computation of Fourier transform on $\mathbb{S}^{2}$ and $\mathrm{SO}(3)$. We refer for further details to Cohen et al. (2018).

5.3 Graph Neural Networks

Graph Neural Networks (GNNs) are the realisation of our Geometric Deep Learning blueprint on graphs leveraging the properties of the permutation group. GNNs are among the most general class of deep learning architectures currently in existence, and as we will see in this text, most other deep learning architectures can be understood as a special case of the GNN with additional geometric structure.

As per our discussion in Section 4.1, we consider a graph to be specified with an adjacency matrix $\mathbf{A}$ and node features $\mathbf{X}$. We will study GNN architectures that are permutation equivariant functions $\mathbf{F}(\mathbf{X},\mathbf{A})$ constructed by applying shared permutation invariant functions $\phi(\mathbf{x}_{u},\mathbf{X}_{\mathcal{N}_{u}})$ over local neighbourhoods. Under various guises, this local function $\phi$ can be referred to as “diffusion”, “propagation”, or “message passing”, and the overall computation of such $\mathbf{F}$ as a “GNN layer”.

The design and study of GNN layers is one of the most active areas of deep learning at the time of writing, making it a landscape that is challenging to navigate. Fortunately, we find that the vast majority of the literature may be derived from only three “flavours” of GNN layers (Figure 17), which we will present here. These flavours govern the extent to which $\phi$ transforms the neighbourhood features, allowing for varying degrees of complexity when modelling interactions across the graph.

In all three flavours, permutation invariance is ensured by aggregating features from $\mathbf{X}_{\mathcal{N}_{u}}$ (potentially transformed, by means of some function $\psi$) with some permutation-invariant function $\bigoplus$, and then updating the features of node $u$, by means of some function $\phi$. Typically,$\psi$ and $\phi$ are learnable, whereas $\bigoplus$ is realised as a nonparametric operation such as sum, mean, or maximum, though it can also be constructed e.g. using recurrent neural networks (Murphy et al., 2018).

Figure 17: A visualisation of the dataflow for the three flavours of GNN layers, $g$. We use the neighbourhood of node $b$ from Figure 10 to illustrate this. Left-to-right: convolutional, where sender node features are multiplied with a constant, $c_{uv}$; attentional, where this multiplier is implicitly computed via an attention mechanism of the receiver over the sender: $\alpha_{uv}=a(\mathbf{x}_{u},\mathbf{x}_{v})$; and message-passing, where vector-based messages are computed based on both the sender and receiver: $\mathbf{m}_{uv}=\psi(\mathbf{x}_{u},\mathbf{x}_{v})$.

In the convolutional flavour (Kipf and Welling, 2016a; Defferrard et al., 2016; Wu et al., 2019), the features of the neighbourhood nodes are directly aggregated with fixed weights,

$$\mathbf{h}_{u}=\phi\left(\mathbf{x}_{u},\bigoplus\limits_{v\in\mathcal{N}_{u}}c_{uv}\psi(\mathbf{x}_{v})\right).$$ (33)

Here, $c_{uv}$ specifies the importance of node $v$ to node $u$’s representation. It is a constant that often directly depends on the entries in ${\bf A}$ representing the structure of the graph. Note that when the aggregation operator $\bigoplus$ is chosen to be the summation, it can be considered as a linear diffusion or position-dependent linear filtering, a generalisation of convolution.In particular, the spectral filters we have seen in Sections 4.4 and 4.6 fall under this category, as they amount to applying fixed local operators (e.g. the Laplacian matrix) to node-wise signals.

In the attentional flavour (Veličković et al., 2018; Monti et al., 2017; Zhang et al., 2018), the interactions are implicit

$$\mathbf{h}_{u}=\phi\left(\mathbf{x}_{u},\bigoplus\limits_{v\in\mathcal{N}_{u}}a(\mathbf{x}_{u},\mathbf{x}_{v})\psi(\mathbf{x}_{v})\right).$$ (34)

Here, $a$ is a learnable self-attention mechanism that computes the importance coefficients $\alpha_{uv}=a(\mathbf{x}_{u},\mathbf{x}_{v})$ implicitly. They are often softmax-normalised across all neighbours. When $\bigoplus$ is the summation, the aggregation is still a linear combination of the neighbourhood node features, but now the weights are feature-dependent.

Finally, the message-passing flavour (Gilmer et al., 2017; Battaglia et al., 2018) amounts to computing arbitrary vectors (“messages”) across edges,

$$\mathbf{h}_{u}=\phi\left(\mathbf{x}_{u},\bigoplus\limits_{v\in\mathcal{N}_{u}}\psi(\mathbf{x}_{u},\mathbf{x}_{v})\right).$$ (35)

Here, $\psi$ is a learnable message function, computing $v$’s vector sent to $u$, and the aggregation can be considered as a form of message passing on the graph.

One important thing to note is a representational containment between these approaches: convolution $\subseteq$ attention $\subseteq$ message-passing. Indeed, attentional GNNs can represent convolutional GNNs by an attention mechanism implemented as a look-up table $a(\mathbf{x}_{u},\mathbf{x}_{v})=c_{uv}$, and both convolutional and attentional GNNs are special cases of message-passing where the messages are only the sender nodes’ features: $\psi(\mathbf{x}_{u},\mathbf{x}_{v})=c_{uv}\psi(\mathbf{x}_{v})$ for convolutional GNNs and $\psi(\mathbf{x}_{u},\mathbf{x}_{v})=a(\mathbf{x}_{u},\mathbf{x}_{v})\psi(\mathbf{x}_{v})$ for attentional GNNs.

This does not imply that message passing GNNs are always the most useful variant; as they have to compute vector-valued messages across edges, they are typically harder to train and require unwieldy amounts of memory. Further, on a wide range of naturally-occurring graphs, the graph’s edges encode for downstream class similarity (i.e. an edge $(u,v)$ implies that $u$ and $v$ are likely to have the same output). For such graphs (often called homophilous), convolutional aggregation across neighbourhoods is often a far better choice, both in terms of regularisation and scalability. Attentional GNNs offer a “middle-ground”: they allow for modelling complex interactions within neighbourhoods while computing only scalar-valued quantities across the edges, making them more scalable than message-passing.

The “three flavour” categorisation presented here is provided with brevity in mind and inevitably neglects a wealth of nuances, insights, generalisations, and historical contexts to GNN models. Importantly, it excludes higher-dimensional GNN based on the Weisfeiler-Lehman hierarchy and spectral GNNs relying on the explicit computation of the graph Fourier transform.

5.4 Deep Sets, Transformers, and Latent Graph Inference

We close the discussion on GNNs by remarking on permutation-equivariant neural network architectures for learning representations of unordered sets. While sets have the least structure among the domains we have discussed in this text, their importance has been recently highlighted by highly-popular architectures such as Transformers (Vaswani et al., 2017) and Deep Sets (Zaheer et al., 2017). In the language of Section 4.1, we assume that we are given a matrix of node features, ${\bf X}$, but without any specified adjacency or ordering information between the nodes. The specific architectures will arise by deciding to what extent to model interactions between the nodes.

Empty edge set

Unordered sets are provided without any additional structure or geometry whatsoever—hence, it could be argued that the most natural way to process them is to treat each set element entirely independently. This translates to a permutation equivariant function over such inputs, which was already introduced in Section 4.1: a shared transformation applied to every node in isolation. Assuming the same notation as when describing GNNs (Section 5.3), such models can be represented as

$$\mathbf{h}_{u}=\psi(\mathbf{x}_{u}),$$

where $\psi$ is a learnable transformation. It may be observed that this is a special case of a convolutional GNN with $\mathcal{N}_{u}=\{u\}$—or, equivalently, ${\bf A}={\bf I}$. Such an architecture is commonly referred to as Deep Sets, in recognition of the work of Zaheer et al. (2017) that have theoretically proved several universal-approximation properties of such architectures. It should be noted that the need to process unordered sets commonly arises in computer vision and graphics when dealing with point clouds; therein, such models are known as PointNets (Qi et al., 2017).

Complete edge set

While assuming an empty edge set is a very efficient construct for building functions over unordered sets, often we would expect that elements of the set exhibit some form of relational structure—i.e., that there exists a latent graph between the nodes. Setting ${\bf A}={\bf I}$ discards any such structure, and may yield suboptimal performance. Conversely, we could assume that, in absence of any other prior knowledge, we cannot upfront exclude any possible links between nodes. In this approach we assume the complete graph, ${\bf A}={\bf 1}{\bf 1}^{\top}$; equivalently, $\mathcal{N}_{u}=\mathcal{V}$. As we do not assume access to any coefficients of interaction, running convolutional-type GNNs over such a graph would amount to:

$$\mathbf{h}_{u}=\phi\left(\mathbf{x}_{u},\bigoplus_{v\in\mathcal{V}}\psi(\mathbf{x}_{v})\right),$$

where the second input, $\bigoplus_{v\in\mathcal{V}}\psi(\mathbf{x}_{v})$ is identical for all nodes $u$, and as such makes the model equivalently expressive to ignoring that input altogether; i.e. the ${\bf A}={\bf I}$ case mentioned above.

This motivates the use of a more expressive GNN flavour, the attentional,

$$\mathbf{h}_{u}=\phi\left(\mathbf{x}_{u},\bigoplus_{v\in\mathcal{V}}a(\mathbf{x}_{u},\mathbf{x}_{v})\psi(\mathbf{x}_{v})\right)$$ (36)

which yields the self-attention operator, the core of the Transformer architecture (Vaswani et al., 2017). Assuming some kind of normalisation over the attentional coefficients (e.g. softmax), we can constrain all the scalars $a(\mathbf{x}_{u},\mathbf{x}_{v})$ to be in the range $[0,1]$; as such, we can think of self-attention as inferring a soft adjacency matrix, $a_{uv}=a(\mathbf{x}_{u},\mathbf{x}_{v})$, as a byproduct of gradient-based optimisation for some downstream task.

The above perspective means that we can pose Transformers exactly as attentional GNNs over a complete graph (Joshi, 2020).However, this is in apparent conflict with Transformers being initially proposed for modelling sequences—the representations of $\mathbf{h}_{u}$ should be mindful of node $u$’s position in the sequence, which complete-graph aggregation would ignore. Transformers address this issue by introducing positional encodings: the node features $\mathbf{x}_{u}$ are augmented to encode node $u$’s position in the sequence, typically as samples from a sine wave whose frequency is dependent on $u$.

On graphs, where no natural ordering of nodes exists, multiple alternatives were proposed to such positional encodings. While we defer discussing these alternatives for later, we note that one promising direction involves a realisation that the positional encodings used in Transformers can be directly related to the discrete Fourier transform (DFT), and hence to the eigenvectors of the graph Laplacian of a “circular grid”. Hence, Transformers’ positional encodings are implicitly representing our assumption that input nodes are connected in a grid. For more general graph structures, one may simply use the Laplacian eigenvectors of the (assumed) graph—an observation exploited by Dwivedi and Bresson (2020) within their empirically powerful Graph Transformer model.

Inferred edge set

Finally, one can try to learn the latent relational structure, leading to some general ${\bf A}$ that is neither ${\bf I}$ nor ${\bf 1}{\bf 1}^{\top}$. The problem of inferring a latent adjacency matrix ${\bf A}$ for a GNN to use (often called latent graph inference) is of high interest for graph representation learning. This is due to the fact that assuming ${\bf A}={\bf I}$ may be representationally inferior, and ${\bf A}={\bf 1}{\bf 1}^{\top}$ may be challenging to implement due to memory requirements and large neighbourhoods to aggregate over. Additionally, it is closest to the “true” problem: inferring an adjacency matrix ${\bf A}$ implies detecting useful structure between the rows of ${\bf X}$, which may then help formulate hypotheses such as causal relations between variables.

Unfortunately, such a framing necessarily induces a step-up in modelling complexity. Specifically, it requires properly balancing a structure learning objective (which is discrete, and hence challenging for gradient-based optimisation) with any downstream task the graph is used for. This makes latent graph inference a highly challenging and intricate problem.

5.5 Equivariant Message Passing Networks

In many applications of Graph Neural Networks, node features (or parts thereof) are not just arbitrary vectors but coordinates of geometric entities. This is the case, for example, when dealing with molecular graphs: the nodes representing atoms may contain information about the atom type as well as its 3D spatial coordinates. It is desirable to process the latter part of the features in a manner that would transform in the same way as the molecule is transformed in space, in other words, be equivariant to the Euclidean group $\mathrm{E}(3)$ of rigid motions (rotations, translations, and reflections) in addition to the standard permutation equivariance discussed before.

To set the stage for our (slightly simplified) analysis, we will make a distinction between node features $\mathbf{f}_{u}\in\mathbb{R}^{d}$ and node spatial coordinates $\mathbf{x}_{u}\in\mathbb{R}^{3}$; the latter are endowed with Euclidean symmetry structure. In this setting, an equivariant layer explicitly transforms these two inputs separately, yielding modified node features $\mathbf{f}^{\prime}_{u}$ and coordinates $\mathbf{x}^{\prime}_{u}$.

We can now state our desirable equivariance property, following the Geometric Deep Learning blueprint. If the spatial component of the input is transformed by ${\mathfrak{g}}\in\mathrm{E}(3)$ (represented as $\rho({\mathfrak{g}})\mathbf{x}=\mathbf{R}\mathbf{x}+\mathbf{b}$, where $\mathbf{R}$ is an orthogonal matrix modeling rotations and reflections, and $\mathbf{b}$ is a translation vector), the spatial component of the output transforms in the same way (as $\mathbf{x}^{\prime}_{u}\mapsto\mathbf{R}\mathbf{x}^{\prime}_{u}+\mathbf{b}$), whereas $\mathbf{f}^{\prime}_{u}$ remains invariant.

Much like the space of permutation equivariant functions we discussed before in the context of general graphs, there exists a vast amount of $\mathrm{E}(3)$-equivariant layers that would satisfy the constraints above—but not all of these layers would be geometrically stable, or easy to implement. In fact, the space of practically useful equivariant layers may well be easily described by a simple categorisation, not unlike our “three flavours” of spatial GNN layers. One elegant solution was suggested by Satorras et al. (2021) in the form of equivariant message passing. Their model operates as follows:

$$\mathbf{f}^{\prime}_{u} = \phi\left(\mathbf{f}_{u},\bigoplus\limits_{v\in\mathcal{N}_{u}}\psi_{\mathrm{f}}(\mathbf{f}_{u},\mathbf{f}_{v},\|\mathbf{x}_{u}-\mathbf{x}_{v}\|^{2})\right),$$

$$\mathbf{x}^{\prime}_{u} = \mathbf{x}_{u}+\sum_{v\neq u}(\mathbf{x}_{u}-\mathbf{x}_{v})\psi_{\mathrm{c}}(\mathbf{f}_{u},\mathbf{f}_{v},\|\mathbf{x}_{u}-\mathbf{x}_{v}\|^{2})$$

where $\psi_{\mathrm{f}}$ and $\psi_{\mathrm{c}}$ are two distinct (learnable) functions. It can be shown that such an aggregation is equivariant under Euclidean transformations of the spatial coordinates. This is due to the fact that the only dependence of $\mathbf{f}^{\prime}_{u}$ on $\mathbf{x}_{u}$ is through the distances $\|\mathbf{x}_{u}-\mathbf{x}_{v}\|^{2}$, and the action of $\mathrm{E}(3)$ necessarily leaves distances between nodes unchanged. Further, the computations of such a layer can be seen as a particular instance of the “message-passing” GNN flavour, hence they are efficient to implement.

To summarise, in contrast to ordinary GNNs, Satorras et al. (2021) enable the correct treatment of ‘coordinates’ for each point in the graph. They are now treated as a member of the $\mathrm{E}(3)$ group, which means the network outputs behave correctly under rotations, reflections and translations of the input.The features, $\mathbf{f}_{u}$, however, are treated in a channel-wise manner and still assumed to be scalars that do not change under these transformations. This limits the type of spatial information that can be captured within such a framework. For example, it may be desirable for some features to be encoded as vectors—e.g. point velocities—which should change direction under such transformations. Satorras et al. (2021) partially alleviate this issue by introducing the concept of velocities in one variant of their architecture. Velocities are a 3D vector property of each point which rotates appropriately. However, this is only a small subspace of the general representations that could be learned with an $\mathrm{E}(3)$ equivariant network. In general, node features may encode tensors of arbitrary dimensionality that would still transform according to $\mathrm{E}(3)$ in a well-defined manner.

Hence, while the architecture discussed above already presents an elegant equivariant solution for many practical input representations, in some cases it may be desirable to explore a broader collection of functions that satisfy the equivariance property. Existing methods dealing with such settings can be categorised into two classes: irreducible representations (of which the previously mentioned layer is a simplified instance) and regular representations. We briefly survey them here, leaving detailed discussion to future work.

Irreducible representations

Irreducible representations build on the finding that all elements of the roto-translation group can be brought into an irreducible form: a vector that is rotated by a block diagonal matrix. Crucially, each of those blocks is a Wigner D-matrix (the aforementioned Fourier basis for Spherical CNNs). Approaches under this umbrella map from one set of irreducible representations to another using equivariant kernels. To find the full set of equivariant mappings, one can then directly solve the equivariance constraint over these kernels. The solutions form a linear combination of equivariant basis matrices derived by Clebsch-Gordan matrices and the spherical harmonics.

Early examples of the irreducible representations approach include Tensor Field Networks (Thomas et al., 2018) and 3D Steerable CNNs (Weiler et al., 2018), both convolutional models operating on point clouds. The $\mathrm{SE}(3)$-Transformer of Fuchs et al. (2020) extends this framework to the graph domain, using an attentional layer rather than convolutional. Further, while our discussion focused on the special case solution of Satorras et al. (2021), we note that the motivation for rotation or translation equivariant predictions over graphs had historically been explored in other fields, including architectures such as Dynamic Graph CNN (Wang et al., 2019b) for point clouds and efficient message passing models for quantum chemistry, such as SchNet (Schütt et al., 2018) and DimeNet (Klicpera et al., 2020).

Regular representations

While the approach of irreducible representations is attractive, it requires directly reasoning about the underlying group representations, which may be tedious, and only applicable to groups that are compact. Regular representation approaches are more general, but come with an additional computational burden – for exact equivariance they require storing copies of latent feature embeddings for all group elements.

One promising approach in this space aims to observe equivariance to Lie groups—through definitions of exponential and logarithmic maps—with the promise of rapid prototyping across various symmetry groups. While Lie groups are out of scope for this section, we refer the reader to two recent successful instances of this direction: the LieConv of Finzi et al. (2020), and the LieTransformer of Hutchinson et al. (2020).

The approaches covered in this section represent popular ways of processing data on geometric graphs in an way that is explicitly equivariant to the underlying geometry. As discussed in Section 4.6, meshes are a special instance of geometric graphs that can be understood as discretisations of continuous surfaces. We will study mesh-specific equivariant neural networks next.

5.6 Intrinsic Mesh CNNs

Meshes, in particular, triangular ones, are the ‘bread and butter’ of computer graphics and perhaps the most common way of modeling 3D objects. The remarkable success of deep learning in general and CNNs in computer vision in particular has lead to a keen interest in the graphics and geometry processing community around the mid-2010sto construct similar architectures for mesh data.

Geodesic patches

Most of the architectures for deep learning on meshes implement convolutional filters of the form (21) by discretising or approximating the exponential map and expressing the filter in a coordinate system of the tangent plane. Shooting a geodesic $\gamma:[0,T]\rightarrow\Omega$ from a point $u=\gamma(0)$ to nearby point $v=\gamma(T)$ defines a local system of geodesic polar coordinates $(r(u,v),\vartheta(u,v))$ where $r$ is the geodesic distance between $u$ and $v$ (length of the geodesic $\gamma$) and $\vartheta$ is the angle between $\gamma^{\prime}(0)$ and some local reference direction. This allows to define a geodesic patch $x(u,r,\vartheta)=x(\exp_{u}\tilde{\omega}(r,\vartheta))$, where $\tilde{\omega}_{u}:[0,R]\times[0,2\pi)\rightarrow T_{u}\Omega$ is the local polar frame.

On a surfacediscretised as a mesh, a geodesic is a poly-line that traverses the triangular faces. Traditionally, geodesics have been computed using the Fast Marching algorithm Kimmel and Sethian (1998), an efficient numerical approximation of a nonlinear PDE called the eikonal equation encountered in physical models of wave propagation in a medium. This scheme was adapted by Kokkinos et al. (2012) for the computation of local geodesic patches and later reused by Masci et al. (2015) for the construction of Geodesic CNNs, the first intrinsic CNN-like architectures on meshes.

Isotropic filters

Importantly, in the definition of the geodesic patch we have ambiguity in the choice of the reference direction and the patch orientation. This is exactly the ambiguity of the choice of the gauge, and our local system of coordinates is defined up to arbitrary rotation (or a shift in the angular coordinate, $x(u,r,\vartheta+\vartheta_{0})$), which can be different at every node. Perhaps the most straightforward solution is to use isotropic filters of the form $\theta(r)$ that perform a direction-independent aggregation of the neighbour features,

$$(x\star\theta)(u)=\int_{0}^{R}\int_{0}^{2\pi}x(u,r,\vartheta)\theta(r)\mathrm{d}r\mathrm{d}\vartheta.$$

Spectral filters discussed in Sections 4.4–4.6 fall under this category: they are based on the Laplacian operator, which is isotropic. Such an approach, however, discards important directional information, and might fail to extract edge-like features.

Fixed gauge

An alternative, to which we have already alluded in Section 4.4, is to fix some gauge. Monti et al. (2017) used the principal curvature directions: while this choice is not intrinsic and may ambiguous at flat points (where curvature vanishes) or uniform curvature (such as on a perfect sphere), the authors showed that it is reasonable for dealing with deformable human body shapes, which are approximately piecewise-rigid. Later works, e.g. Melzi et al. (2019), showed reliable intrinsic construction of gauges on meshes, computed as (intrinsic) gradients of intrinsic functions. While such tangent fields might have singularities (i.e., vanish at some points), the overall procedure is very robust to noise and remeshing.

Angular pooling

Another approach, referred to as angular max pooling, was used by Masci et al. (2015). In this case, the filter $\theta(r,\vartheta)$ is anisotropic, but its matching with the function is performed over all the possible rotations, which are then aggregated:

$$(x\star\theta)(u)=\max_{\vartheta_{0}\in[0,2\pi)}\,\,\int_{0}^{R}\int_{0}^{2\pi}x(u,r,\vartheta)\theta(r,\vartheta+\vartheta_{0})\mathrm{d}r\mathrm{d}\vartheta.$$

Conceptually, this can be visualised as correlating geodesic patches with a rotating filter and collecting the strongest responses.

On meshes, the continuous integrals can be discretised using a construction referred to as patch operators (Masci et al., 2015). In a geodesic patch around node $u$, the neighbour nodes $\mathcal{N}_{u}$,represented in the local polar coordinates as $(r_{uv},\vartheta_{uv})$, are weighted by a set of weighting functions $w_{1}(r,\vartheta),\ldots,w_{K}(r,\vartheta)$ (shown in Figure 18 and acting as ‘soft pixels’) and aggregated,

$$(x\star\theta)_{u}=\frac{\sum_{k=1}^{K}w_{k}\sum_{v\in\mathcal{N}_{u}}(r_{uv},\vartheta_{uv})x_{v}\,\theta_{k}}{\sum_{k=1}^{K}w_{k}\sum_{v\in\mathcal{N}_{u}}(r_{uv},\vartheta_{uv})\theta_{k}}$$

(here $\theta_{1},\ldots,\theta_{K}$ are the learnable coefficients of the filter). Multi-channel features are treated channel-wise, with a family of appropriate filters. Masci et al. (2015); Boscaini et al. (2016a) used pre-defined weighting functions $w$, while Monti et al. (2017) further allowed them to be learnable.

Figure 18: Left-to-right: examples of patch operators used in Geodesic CNN (Masci et al., 2015), Anisotropic CNN (Boscaini et al., 2016b) and MoNet (Monti et al., 2017), with the level sets of the weighting functions $w_{k}(r,\vartheta)$ shown in red.

Gauge-equivariant filters

Both isotropic filters and angular max pooling lead to features that are invariant to gauge transformations; they transform according to the trivial representation $\rho({\mathfrak{g}})=1$ (where ${\mathfrak{g}}\in\mathrm{SO}(2)$ is a rotation of the local coordinate frame). This point of view suggests another approach, proposed by Cohen et al. (2019); de Haan et al. (2020) and discussed in Section 4.5, where the features computed by the network are associated with an arbitrary representation $\rho$ of the structure group ${\mathfrak{G}}$ (e.g. $\mathrm{SO(2)}$ or $\mathrm{O(2)}$ of rotations or rotations+reflections of the coordinate frame, respectively). Tangent vectors transform according to the standard representation $\rho({\mathfrak{g}})={\mathfrak{g}}$. As another example, the feature vector obtained by matching $n$ rotated copies of the same filter transforms by cyclic shifts under rotations of the gauge; this is known as the regular representation of the cyclic group $C_{n}$.

As discussed in Section 4.5, when dealing with such geometric features (associated to a non-trivial representation), we must first parallel transport them to the same vector space before applying the filter. On a mesh, this can be implemented via the following message passing mechanism described by de Haan et al. (2020). Let $\mathbf{x}_{u}\in\mathbb{R}^{d}$ be a $d$-dimensional input feature at mesh node $u$. This feature is expressed relative to an (arbitrary) choice of gauge at $u$, and is assumed to transform according to a representation $\rho_{\textup{in}}$ of ${\mathfrak{G}}=\operatorname{SO}(2)$ under rotations of the gauge. Similarly, the output features $\mathbf{h}_{u}$ of the mesh convolution are $d^{\prime}$ dimensional and should transform according to $\rho_{\textup{out}}$ (which can be chosen at will by the network designer).

By analogy to Graph Neural Networks, we can implement the gauge-equivariant convolution (23) on meshes by sending messages from the neighbours $\mathcal{N}_{u}$ of $u$ (and from $u$ itself) to $u$:

$$\mathbf{h}_{u}=\bm{\Theta}_{\textup{self}}\;\mathbf{x}_{u}+\sum_{v\in\mathcal{N}_{u}}\bm{\Theta}_{\textup{neigh}}(\vartheta_{uv})\rho({\mathfrak{g}}_{v\rightarrow u})\mathbf{x}_{v},$$ (37)

where $\bm{\Theta}_{\textup{self}},\bm{\Theta}_{\textup{neigh}}(\vartheta_{uv})\in\mathbb{R}^{d^{\prime}\times d}$ are learned filter matrices. The structure group element ${\mathfrak{g}}_{v\rightarrow u}\in\operatorname{SO}(2)$ denotes the effect of parallel transport from $v$ to $u$, expressed relative to the gauges at $u$ and $v$, and can be precomputed for each mesh. Its action is encoded by a transporter matrix $\rho({\mathfrak{g}}_{v\rightarrow u})\in\mathbb{R}^{d\times d}$.The matrix $\bm{\Theta}_{\textup{neigh}}(\vartheta_{uv})$ depends on the angle $\vartheta_{uv}$ of the neighbour $v$ to the reference direction (e.g. first axis of the frame) at $u$, so this kernel is anisotropic: different neighbours are treated differently.

As explained in Section 4.5, for $\mathbf{h}(u)$ to be a well-defined geometric quantity, it should transform as $\mathbf{h}(u)\mapsto\rho_{\textup{out}}({\mathfrak{g}}^{-1}(u))\mathbf{h}(u)$ under gauge transformations. This will be the case when $\bm{\Theta}_{\textup{self}}\rho_{\textup{in}}(\vartheta)=\rho_{\textup{out}}(\vartheta)\bm{\Theta}_{\textup{self}}$ for all $\vartheta\in\operatorname{SO}(2)$,and $\bm{\Theta}_{\textup{neigh}}(\vartheta_{uv}-\vartheta)\rho_{\textup{in}}(\vartheta)=\rho_{\textup{out}}(\vartheta)\bm{\Theta}_{\textup{neigh}}(\vartheta_{uv})$. Since these constraints are linear, the space of matrices $\bm{\Theta}_{\textup{self}}$ and matrix-valued functions $\bm{\Theta}_{\textup{neigh}}$ satisfying these constraints is a linear subspace, and so we can parameterise them as a linear combination of basis kernels with learnable coefficients: $\bm{\Theta}_{\textup{self}}=\sum_{i}\alpha_{i}\bm{\Theta}_{\textup{self}}^{i}$ and $\bm{\Theta}_{\textup{neigh}}=\sum_{i}\beta_{i}\bm{\Theta}_{\textup{neigh}}^{i}$.

5.7 Recurrent Neural Networks

Our discussion has thus far always assumed the inputs to be solely spatial across a given domain. However, in many common use cases, the inputs can also be considered sequential (e.g. video, text or speech). In this case, we assume that the input consists of arbitrarily many steps, wherein at each step $t$ we are provided with an input signal, which we represent as $\mathbf{X}^{(t)}\in\mathcal{X}(\Omega^{(t)})$.

While in general the domain can evolve in time together with the signals on it, it is typically assumed that the domain is kept fixed across all the $t$, i.e. $\Omega^{(t)}=\Omega$. Here, we will exclusively focus on this case, but note that exceptions are common. Social networks are an example where one often has to account for the domain changing through time, as new links are regularly created as well as erased. The domain in this setting is often referred to as a dynamic graph (Xu et al., 2020a; Rossi et al., 2020).

Often, the individual $\mathbf{X}^{(t)}$ inputs will exhibit useful symmetries and hence may be nontrivially treated by any of our previously discussed architectures. Some common examples include: videos ($\Omega$ is a fixed grid, and signals are a sequence of frames); fMRI scans ($\Omega$ is a fixed mesh representing the geometry of the brain cortex, where different regions are activated at different times as a response to presented stimuli); and traffic flow networks ($\Omega$ is a fixed graph representing the road network, on which e.g. the average traffic speed is recorded at various nodes).

Let us assume an encoder function $f(\mathbf{X}^{(t)})$ providing latent representations at the level of granularity appropriate for the problem and respectful of the symmetries of the input domain. As an example, consider processing video frames: that is, at each timestep, we are given a grid-structured input represented as an $n\times d$ matrix $\mathbf{X}^{(t)}$, where $n$ is the number of pixels (fixed in time) and $d$ is the number of input channels (e.g. $d=3$ for RGB frames). Further, we are interested in analysis at the level of entire frames, in which case it is appropriate to implement $f$ as a translation invariant CNN, outputting a $k$-dimensional representation $\mathbf{z}^{(t)}=f(\mathbf{X}^{(t)})$ of the frame at time-step $t$.

We are now left with the task of appropriately summarising a sequence of vectors $\mathbf{z}^{(t)}$ across all the steps. A canonical way to dynamically aggregate this information in a way that respects the temporal progression of inputs and also easily allows for online arrival of novel data-points, is using a Recurrent Neural Network (RNN). What we will show here is that RNNs are an interesting geometric architecture to study in their own right, since they implement a rather unusual type of symmetry over the inputs $\mathbf{z}^{(t)}$.

SimpleRNNs

Figure 19: Illustration of processing video input with RNNs. Each input video frame $\mathbf{X}^{(t)}$ is processed using a shared function $f$—e.g. a translation invariant CNN—into a flat representation $\mathbf{z}^{(t)}$. Then the RNN update function $R$ is iterated across these vectors, iteratively updating a summary vector $\mathbf{h}^{(t)}$ which summarises all the inputs up to and including $\mathbf{z}^{(t)}$. The computation is seeded with an initial summary vector $\mathbf{h}^{(0)}$, which may be either pre-determined or learnable.

At each step, the recurrent neural network computes an $m$-dimensional summary vector $\mathbf{h}^{(t)}$ of all the input steps up to and including $t$. This (partial) summary is computed conditional on the current step’s features and the previous step’s summary, through a shared update function, $R:\mathbb{R}^{k}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$, as follows (see Figure 19 for a summary):

$$\mathbf{h}^{(t)}=R(\mathbf{z}^{(t)},\mathbf{h}^{(t-1)})$$ (38)

and, as both $\mathbf{z}^{(t)}$ and $\mathbf{h}^{(t-1)}$ are flat vector representations, $R$ may be most easily expressed as a single fully-connected neural network layer (often known as SimpleRNN; see Elman (1990); Jordan (1997)):

$$\mathbf{h}^{(t)}=\sigma(\mathbf{W}\mathbf{z}^{(t)}+\mathbf{U}\mathbf{h}^{(t-1)}+\mathbf{b})$$ (39)

where $\mathbf{W}\in\mathbb{R}^{k\times m}$, $\mathbf{U}\in\mathbb{R}^{m\times m}$ and $\mathbf{b}\in\mathbb{R}^{m}$ are learnable parameters, and $\sigma$ is an activation function. While this introduces loops in the network’s computational graph, in practice the network is unrolled for an appropriate number of steps, allowing for backpropagation through time (Robinson and Fallside, 1987; Werbos, 1988; Mozer, 1989) to be applied.

The summary vectors may then be appropriately leveraged for the downstream task—if a prediction is required at every step of the sequence, then a shared predictor may be applied to each $\mathbf{h}^{(t)}$ individually. For classifying entire sequences, typically the final summary, $\mathbf{h}^{(T)}$, is passed to a classifier. Here, $T$ is the length of the sequence.

Specially, the initial summary vector is usually either set to the zero-vector, i.e. $\mathbf{h}^{(0)}=\mathbf{0}$, or it is made learnable. Analysing the manner in which the initial summary vector is set also allows us to deduce an interesting form of translation equivariance exhibited by RNNs.

Translation equivariance in RNNs

Since we interpret the individual steps $t$ as discrete time-steps, the input vectors $\mathbf{z}^{(t)}$ can be seen as living on a one-dimensionalgrid of time-steps. While it might be attractive to attempt extending our translation equivariance analysis from CNNs here, it cannot be done in a trivial manner.

To see why, let us assume that we have produced a new sequence $\mathbf{z}^{\prime(t)}=\mathbf{z}^{(t+1)}$ by performing a left-shift of our sequence by one step. It might be tempting to attempt showing $\mathbf{h}^{\prime(t)}=\mathbf{h}^{(t+1)}$, as one expects with translation equivariance; however, this will not generally hold. Consider $t=1$; directly applying and expanding the update function, we recover the following:

$$\mathbf{h}^{\prime(1)} =R(\mathbf{z}^{\prime(1)},\mathbf{h}^{(0)})=R(\mathbf{z}^{(2)},\mathbf{h}^{(0)})$$ (40)

$$\mathbf{h}^{(2)} =R(\mathbf{z}^{(2)},\mathbf{h}^{(1)})=R(\mathbf{z}^{(2)},R(\mathbf{z}^{(1)},\mathbf{h}^{(0)}))$$ (41)

Hence, unless we can guarantee that $\mathbf{h}^{(0)}=R(\mathbf{z}^{(1)},\mathbf{h}^{(0)})$, we will not recover translation equivariance. Similar analysis can then be done for steps $t>1$.

Fortunately, with a slight refactoring of how we represent $\mathbf{z}$, and for a suitable choice of $R$, it is possible to satisfy the equality above, and hence demonstrate a setting in which RNNs are equivariant to shifts. Our problem was largely one of boundary conditions: the equality above includes $\mathbf{z}^{(1)}$, which our left-shift operation destroyed. To abstract this problem away, we will observe how an RNN processes an appropriately left-padded sequence, $\bar{\mathbf{z}}^{(t)}$, defined as follows:

$$\bar{\mathbf{z}}^{(t)}=\begin{cases}\mathbf{0}&t\leq t^{\prime}\\ \mathbf{z}^{(t-t^{\prime})}&t>t^{\prime}\end{cases}$$

Such a sequence now allows for left-shiftingby up to $t^{\prime}$ steps without destroying any of the original input elements. Further, note we do not need to handle right-shifting separately; indeed, equivariance to right shifts naturally follows from the RNN equations.

We can now again analyse the operation of the RNN over a left-shifted verson of $\bar{\mathbf{z}}^{(t)}$, which we denote by $\bar{\mathbf{z}}^{\prime(t)}=\bar{\mathbf{z}}^{(t+1)}$, as we did in Equations 40–41:

$$\mathbf{h}^{\prime(1)} =R(\bar{\mathbf{z}}^{\prime(1)},\mathbf{h}^{(0)})=R(\bar{\mathbf{z}}^{(2)},\mathbf{h}^{(0)})$$

$$\mathbf{h}^{(2)} =R(\bar{\mathbf{z}}^{(2)},\mathbf{h}^{(1)})=R(\bar{\mathbf{z}}^{(2)},R(\bar{\mathbf{z}}^{(1)},\mathbf{h}^{(0)}))=R(\bar{\mathbf{z}}^{(2)},R(\mathbf{0},\mathbf{h}^{(0)}))$$

where the substitution $\bar{\mathbf{z}}^{(1)}=\mathbf{0}$ holds as long as $t^{\prime}\geq 1$, i.e. as long as any padding is applied. Now, we can guarantee equivariance to left-shifting by one step ($\mathbf{h}^{\prime(t)}=\mathbf{h}^{(t+1)}$) as long as $\mathbf{h}^{(0)}=R(\mathbf{0},\mathbf{h}^{(0)})$.

Said differently, $\mathbf{h}^{(0)}$ must be chosen to be a fixed point of a function $\gamma(\mathbf{h})=R(\mathbf{0},\mathbf{h})$. If the update function $R$ is conveniently chosen, then not only can we guarantee existence of such fixed points, but we can even directly obtain them by iterating the application of $R$ until convergence; e.g., as follows:

$$\mathbf{h}_{0}=\mathbf{0}\qquad\mathbf{h}_{k+1}=\gamma(\mathbf{h}_{k}),$$ (42)

where the index $k$ refers to the iteration of $R$ in our computation, as opposed to the the index $(t)$ denoting the time step of the RNN. If we choose $R$ such that $\gamma$ is a contraction mapping, such an iteration will indeed converge to a unique fixed point. Accordingly, we can then iterate Equation (42) until $\mathbf{h}_{k+1}=\mathbf{h}_{k}$, and we can set $\mathbf{h}^{(0)}=\mathbf{h}_{k}$. Note that this computation is equivalent to left-padding the sequence with “sufficiently many” zero-vectors.

Depth in RNNs

It is also easy to stack multiple RNNs—simply use the $\mathbf{h}^{(t)}$ vectors as an input sequence for a second RNN. This kind of construction is occasionally called a “deep RNN”, which is potentially misleading. Effectively, due to the repeated application of the recurrent operation, even a single RNN “layer” has depth equal to the number of input steps.

This often introduces uniquely challenging learning dynamics when optimising RNNs, as each training example induces many gradient updates to the shared parameters of the update network. Here we will focus on perhaps the most prominent such issue—that of vanishing and exploding gradients (Bengio et al., 1994)—which is especially problematic in RNNs, given their depth and parameter sharing. Further, it has single-handedly spurred some of the most influential research on RNNs. For a more detailed overview, we refer the reader to Pascanu et al. (2013), who have studied the training dynamics of RNNs in great detail, and exposed these challenges from a variety of perspectives: analytical, geometrical, and the lens of dynamical systems.

To illustrate vanishing gradients, consider a SimpleRNN with a sigmoidal activation function $\sigma$, whose derivative magnitude $|\sigma^{\prime}|$ is always between $0$ and $1$. Multiplying many such values results in gradients that quickly tend to zero, implying that early steps in the input sequence may not be able to have influence in updating the network parameters at all.

For example, consider the next-word prediction task (common in e.g. predictive keyboards), and the input text “Petar is Serbian. He was born on …[long paragraph] …Petar currently lives in  ”. Here, predicting the next word as “Serbia” may only be reasonably concluded by considering the very start of the paragraph—but gradients have likely vanished by the time they reach this input step, making learning from such examples very challenging.

Deep feedforward neural networks have also suffered from the vanishing gradient problem, until the invention of the ReLU activation (which has gradients equal to exactly zero or one—thus fixing the vanishing gradient problem). However, in RNNs, using ReLUs may easily lead to exploding gradients, as the output space of the update function is now unbounded, and gradient descent will update the cell once for every input step, quickly building up the scale of the updates. Historically, the vanishing gradient phenomenon was recognised early on as a significant obstacle in the use of recurrent networks. Coping with this problem motivated the development of more sophisticated RNN layers, which we describe next.

5.8 Long Short-Term Memory networks

Figure 20: The dataflow of the long short-term memory (LSTM), with its components and memory cell ($M$) clearly highlighted. Based on the current input $\mathbf{z}^{(t)}$, previous summary $\mathbf{h}^{(t-1)}$ and previous cell state $\mathbf{c}^{(t-1)}$, the LSTM predicts the updated cell state $\mathbf{c}^{(t)}$ and summary $\mathbf{h}^{(t)}$.

A key invention that significantly reduced the effects of vanishing gradients in RNNs is that of gating mechanisms, which allow the network to selectively overwrite information in a data-driven way. Prominent examples of these gated RNNs include the Long Short-Term Memory (LSTM; Hochreiter and Schmidhuber (1997)) and the Gated Recurrent Unit (GRU; Cho et al. (2014)). Here we will primarily discuss the LSTM—specifically, the variant presented by Graves (2013)—in order to illustrate the operations of such models. Concepts from LSTMs easily carry over to other gated RNNs.

Throughout this section, it will likely be useful to refer to Figure 20, which illustrates all of the LSTM operations that we will discuss in text.

The LSTM augments the recurrent computation by introducing a memory cell, which stores cell state vectors, $\mathbf{c}^{(t)}\in\mathbb{R}^{m}$, that are preserved between computational steps. The LSTM computes summary vectors, $\mathbf{h}^{(t)}$, directly based on $\mathbf{c}^{(t)}$, and $\mathbf{c}^{(t)}$ is, in turn, computed using $\mathbf{z}^{(t)}$, $\mathbf{h}^{(t-1)}$ and $\mathbf{c}^{(t-1)}$. Critically, the cell is not completely overwritten based on $\mathbf{z}^{(t)}$ and $\mathbf{h}^{(t-1)}$, which would expose the network to the same issues as the SimpleRNN. Instead, a certain quantity of the previous cell state may be retained—and the proportion by which this occurs is explicitly learned from data.

Just like in SimpleRNN, we compute features by using a single fully-connected neural network layer over the current input step and previous summary:

$$\widetilde{\mathbf{c}}^{(t)}=\tanh(\mathbf{W}_{c}\mathbf{z}^{(t)}+\mathbf{U}_{c}\mathbf{h}^{(t-1)}+\mathbf{b}_{c})$$ (43)

But, as mentioned, we do not allow all of this vector to enter the cell—hence why we call it the vector of candidate features, and denote it as $\widetilde{\mathbf{c}}^{(t)}$. Instead, the LSTM directly learns gating vectors, which are real-valued vectors in the range $[0,1]$, and decide how much of the signal should be allowed to enter, exit, and overwrite the memory cell.

Three such gates are computed, all based on $\mathbf{z}^{(t)}$ and $\mathbf{h}^{(t-1)}$: the input gate $\mathbf{i}^{(t)}$, which computes the proportion of the candidate vector allowed to enter the cell; the forget gate $\mathbf{f}^{(t)}$, which computes the proportion of the previous cell state to be retained, and the output gate $\mathbf{o}^{(t)}$, which computes the proportion of the new cell state to be used for the final summary vector. Typically all of these gates are also derived using a single fully connected layer, albeit with the logistic sigmoid activation $\mathrm{logistic}(x)=\frac{1}{1+\exp(-x)}$, in order to guarantee that the outputs are in the $[0,1]$ range:

$$\mathbf{i}^{(t)} =\mathrm{logistic}(\mathbf{W}_{\mathrm{i}}\mathbf{z}^{(t)}+\mathbf{U}_{\mathrm{i}}\mathbf{h}^{(t-1)}+\mathbf{b}_{\mathrm{i}})$$ (44)

$$\mathbf{f}^{(t)} =\mathrm{logistic}(\mathbf{W}_{\mathrm{f}}\mathbf{z}^{(t)}+\mathbf{U}_{\mathrm{f}}\mathbf{h}^{(t-1)}+\mathbf{b}_{\mathrm{f}})$$ (45)

$$\mathbf{o}^{(t)} =\mathrm{logistic}(\mathbf{W}_{\mathrm{o}}\mathbf{z}^{(t)}+\mathbf{U}_{\mathrm{o}}\mathbf{h}^{(t-1)}+\mathbf{b}_{\mathrm{o}})$$ (46)

Finally, these gates are appropriately applied to decode the new cell state, $\mathbf{c}^{(t)}$, which is then modulated by the output gate to produce the summary vector $\mathbf{h}^{(t)}$, as follows:

$$\mathbf{c}^{(t)} =\mathbf{i}^{(t)}\odot\widetilde{\mathbf{c}}^{(t)}+\mathbf{f}^{(t)}\odot\mathbf{c}^{(t-1)}$$ (47)

$$\mathbf{h}^{(t)} =\mathbf{o}^{(t)}\odot\tanh(\mathbf{c}^{(t)})$$ (48)

where $\odot$ is element-wise vector multiplication. Applied together, Equations (43)–(48) completely specify the update rule for the LSTM, which now takes into account the cell vector $\mathbf{c}^{(t)}$ as well:

$$(\mathbf{h}^{(t)},\mathbf{c}^{(t)})=R(\mathbf{z}^{(t)},(\mathbf{h}^{(t-1)},\mathbf{c}^{(t-1)}))$$

Note that, as the values of $\mathbf{f}^{(t)}$ are derived from $\mathbf{z}^{(t)}$ and $\mathbf{h}^{(t-1)}$—and therefore directly learnable from data—the LSTM effectively learns how to appropriately forget past experiences. Indeed, the values of $\mathbf{f}^{(t)}$ directly appear in the backpropagation update for all the LSTM parameters ($\mathbf{W}_{*},\mathbf{U}_{*},\mathbf{b}_{*}$), allowing the network to explicitly control, in a data-driven way, the degree of vanishing for the gradients across the time steps.

Besides tackling the vanishing gradient issue head-on, it turns out that gated RNNs also unlock a very useful form of invariance to time-warping transformations, which remains out of reach of SimpleRNNs.

Time warping invariance of gated RNNs

We will start by illustrating, in a continuous-time setting, what does it mean to warp time, and what is required of a recurrent model in order to achieve invariance to such transformations. Our exposition will largely follow the work of Tallec and Ollivier (2018), that initially described this phenomenon—and indeed, they were among the first to actually study RNNs from the lens of invariances.

Let us assume a continuous time-domain signal $z(t)$, on which we would like to apply an RNN. To align the RNN’s discrete-time computation of summary vectors $\mathbf{h}^{(t)}$with an analogue in the continuous domain, $h(t)$, we will observe its linear Taylor expansion:

$$h(t+\delta)\approx h(t)+\delta\frac{\mathrm{d}h(t)}{\mathrm{d}t}$$ (49)

and, setting $\delta=1$, we recover a relationship between $h(t)$ and $h(t+1)$, which is exactly what the RNN update function $R$ (Equation 38) computes. Namely, the RNN update function satisfies the following differential equation:

$$\frac{\mathrm{d}h(t)}{\mathrm{d}t}=h(t+1)-h(t)=R(z(t+1),h(t))-h(t)$$ (50)

We would like the RNN to be resilient to the way in which the signal is sampled (e.g. by changing the time unit of measurement), in order to account for any imperfections or irregularities therein. Formally, we denote a time warpingoperation $\tau:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$, as any monotonically increasing differentiable mapping between times. The notation $\tau$ is chosen because time warping represents an automorphism of time.

Further, we state that a class of models is invariant to time warping if, for any model of the class and any such $\tau$, there exists another (possibly the same) model from the class that processes the warped data in the same way as the original model did in the non-warped case.

This is a potentially very useful property. If we have an RNN class capable of modelling short-term dependencies well, and we can also show that this class is invariant to time warping, then we know it is possible to train such a model in a way that will usefully capture long-term dependencies as well (as they would correspond to a time dilation warping of a signal with short-term dependencies). As we will shortly see, it is no coincidence that gated RNN models such as the LSTM were proposed to model long-range dependencies. Achieving time warping invariance is tightly coupled with presence of gating mechanisms, such as the input/forget/output gates of LSTMs.

When time gets warped by $\tau$, the signal observed by the RNN at time $t$ is $z(\tau(t))$ and, to remain invariant to such warpings, it should predict an equivalently-warped summary function $h(\tau(t))$. Using Taylor expansion arguments once more, we derive a form of Equation 50 for the warped time, that the RNN update $R$ should satisfy:

$$\frac{\mathrm{d}h(\tau(t))}{\mathrm{d}\tau(t)}=R(z(\tau(t+1)),h(\tau(t)))-h(\tau(t))$$ (51)

However, the above derivative is computed with respect to the warped time $\tau(t)$, and hence does not take into account the original signal. To make our model take into account the warping transformation explicitly, we need to differentiate the warped summary function with respect to $t$. Applying the chain rule, this yields the following differential equation:

$$\frac{\mathrm{d}h(\tau(t))}{\mathrm{d}t}=\frac{\mathrm{d}h(\tau(t))}{\mathrm{d}\tau(t)}\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}=\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}R(z(\tau(t+1)),h(\tau(t)))-\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}h(\tau(t))$$ (52)

and, for our (continuous-time) RNN to remain invariant to any time warping $\tau(t)$, it needs to be able to explicitly represent the derivative $\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}$, which is not assumed known upfront! We need to introduce a learnable function $\Gamma$ which approximates this derivative. For example, $\Gamma$ could be a neural network taking into account $z(t+1)$ and $h(t)$ and predicting scalar outputs.

Now, remark that, from the point of view of a discrete RNN model under time warping, its input $\mathbf{z}^{(t)}$ will correspond to $z(\tau(t))$, and its summary $\mathbf{h}^{(t)}$ will correspond to $h(\tau(t))$. To obtain the required relationship of $\mathbf{h}^{(t)}$ to $\mathbf{h}^{(t+1)}$ in order to remain invariant to time warping, we will use a one-step Taylor expansion of $h(\tau(t))$:

$$h(\tau(t+\delta))\approx h(\tau(t))+\delta\frac{\mathrm{d}h(\tau(t))}{\mathrm{d}t}$$

and, once again, setting $\delta=1$ and substituting Equation 52, then discretising:

$$\mathbf{h}^{(t+1)} =\mathbf{h}^{(t)}+\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}R(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})-\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}\mathbf{h}^{(t)}$$

$$=\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}R(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})+\left(1-\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}\right)\mathbf{h}^{(t)}$$

Finally, we swap $\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}$ with the aforementioned learnable function, $\Gamma$. This gives us the required form for our time warping-invariant RNN:

$$\mathbf{h}^{(t+1)}=\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})R(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})+(1-\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)}))\mathbf{h}^{(t)}$$ (53)

We may quickly deduce that SimpleRNNs (Equation 39) are not time warping invariant, given that they do not feature the second term in Equation 53. Instead, they fully overwrite $\mathbf{h}^{(t)}$ with $R(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})$, which corresponds to assuming no time warping at all; $\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}=1$, i.e. $\tau(t)=t$.

Further, our link between continuous-time RNNs and the discrete RNN based on $R$ rested on the accuracy of the Taylor approximation, which holds only if the time-warping derivative is not too large, i.e., $\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}\lesssim 1$. The intuitive explanation of this is: if our time warping operation ever contracts time in a way that makes time increments ($t\rightarrow t+1$) large enough that intermediate data changes are not sampled, the model can never hope to process time-warped inputs in the same way as original ones—it simply would not have access to the same information. Conversely, time dilations of any form (which, in discrete terms, correspond to interspersing the input time-series with zeroes) are perfectly allowed within our framework.

Combined with our requirement of monotonically increasing $\tau$ ($\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}>0$), we can bound the output space of $\Gamma$ as $0<\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})<1$, which motivates the use of the logistic sigmoid activation for $\Gamma$, e.g.:

$$\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})=\mathrm{logistic}({\bf W}_{\Gamma}\mathbf{z}^{(t+1)}+{\bf U}_{\Gamma}\mathbf{h}^{(t)}+\mathbf{b}_{\Gamma})$$

exactly matching the LSTM gating equations (e.g. Equation 44). The main difference is that LSTMs compute gating vectors, whereas Equation 53 implies $\Gamma$ should output a scalar. Vectorised gates (Hochreiter, 1991) allow to fit a different warping derivative in every dimension of $\mathbf{h}^{(t)}$, allowing for reasoning over multiple time horizons simultaneously.

It is worth taking a pause here to summarise what we have done. By requiring that our RNN class is invariant to (non-destructive) time warping, we have derived the necessary form that it must have (Equation 53), and showed that it exactly corresponds to the class of gated RNNs. The gates’ primary role under this perspective is to accurately fit the derivative $\frac{\mathrm{d}\tau(t)}{\mathrm{d}t}$ of the warping transformation.

The notion of class invariance is somewhat distinct from the invariances we studied previously. Namely, once we train a gated RNN on a time-warped input with $\tau_{1}(t)$, we typically cannot zero-shot transferit to a signal warped by a different $\tau_{2}(t)$. Rather, class invariance only guarantees that gated RNNs are powerful enough to fit both of these signals in the same manner, but potentially with vastly different model parameters. That being said, the realisation that effective gating mechanisms are tightly related to fitting the warping derivative can yield useful prescriptions for gated RNN optimisation, as we now briefly demonstrate.

For example, we can often assume that the range of the dependencies we are interested in tracking within our signal will be in the range $[T_{l},T_{h}]$ time-steps.

By analysing the analytic solutions to Equation 52, it can be shown that the characteristic forgetting time of $\mathbf{h}^{(t)}$ by our gated RNN is proportional to $\frac{1}{\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})}$. Hence, we would like our gating values to lie between $\left[\frac{1}{T_{h}},\frac{1}{T_{m}}\right]$ in order to effectively remember information within the assumed range.

Further, if we assume that $\mathbf{z}^{(t)}$ and $\mathbf{h}^{(t)}$ are roughly zero-centered—which is a common by-product of applying transformations such as layer normalisation (Ba et al., 2016)—we can assume that $\mathbb{E}[\Gamma(\mathbf{z}^{(t+1)},\mathbf{h}^{(t)})]\approx\mathrm{logistic}(\mathbf{b}_{\Gamma})$. Controlling the bias vector of the gating mechanism is hence a very powerful way of controlling the effective gate value.

Combining the two observations, we conclude that an appropriate range of gating values can be obtained by initialising $\mathbf{b}_{\Gamma}\sim-\log(\mathcal{U}(T_{l},T_{h})-1)$, where $\mathcal{U}$ is the uniform real distribution. Such a recommendation was dubbed chrono initialisation by Tallec and Ollivier (2018), and has been empirically shown to improve the long-range dependency modelling of gated RNNs.

Sequence-to-sequence learning with RNNs

One prominent historical example of using RNN-backed computation are sequence-to-sequence translation tasks, such as machine translation of natural languages. The pioneering seq2seq work by Sutskever et al. (2014) achieved this by passing the summary vector, $\mathbf{h}^{(T)}$ as an initial input for a decoder RNN, with outputs of RNN blocks being given as inputs for the next step.

Figure 21: One typical example of a seq2seq architecture with an RNN encoder $R_{\mathrm{enc}}$ and RNN decoder $R_{\mathrm{dec}}$. The decoder is seeded with the final summary vector $\mathbf{h}^{(T)}$ coming out of the encoder, and then proceeds in an autoregressive fashion: at each step, the predicted output from the previous step is fed back as input to $R_{\mathrm{dec}}$. The bottleneck problem is also illustrated with the red lines: the summary vector $\mathbf{h}^{(T)}$ is pressured to store all relevant information for translating the input sequence, which becomes increasingly challenging as the input length grows.

This placed substantial representational pressure on the summary vector, $\mathbf{h}^{(T)}$. Within the context of deep learning, $\mathbf{h}^{(T)}$ is sometimes referred to as a bottleneck. Its fixed capacity must be sufficient for representing the content of the entire input sequence, in a manner that is conducive to generating a corresponding sequence, while also supporting input sequences of substantially different lengths (Figure 21).

In reality, different steps of the output may wish to focus (attend) on different parts of the input, and all such choices are difficult to represent via a bottleneck vector. Following from this observation, the popular recurrent attention model was proposed by Bahdanau et al. (2014). At every step of processing, a query vector is generated by an RNN; this query vector then interacts with the representation of every time-step $\mathbf{h}^{(t)}$, primarily by computing a weighted sum over them. This model pioneered neural content-based attention and predates the success of the Transformer model.

Lastly, while attending offers a soft way to dynamically focus on parts of the input content, substantial work also learnt more explicit ways to direct attention to the input. A powerful algorithmically grounded way of doing so is the pointer network of Vinyals et al. (2015), which proposes a simple modification of recurrent attention to allow for pointing over elements of variable-sized inputs. These findings have then been generalised to the set2set architecture (Vinyals et al., 2016), which generalises seq2seq models to unordered sets, supported by pointer network-backed LSTMs.

6 Problems and Applications

Invariances and symmetries arise all too commonly across data originating in the real world. Hence, it should come as no surprise that some of the most popular applications of machine learning in the 21st century have come about as a direct byproduct of Geometric Deep Learning, perhaps sometimes without fully realising this fact. We would like to provide readers with an overview—by no means comprehensive—of influential works in Geometric Deep Learning and exciting and promising new applications. Our motivation is twofold: to demonstrate specific instances of scientific and industrial problems where the five geometric domains commonly arise, and to serve additional motivation for further study of Geometric Deep Learning principles and architectures.

Chemistry and Drug Design

One of the most promising applications of representation learning on graphs is in computational chemistry and drug development.Traditional drugs are small molecules that are designed to chemically attach (‘bind’) to some target molecule, typically a protein, in order to activate or disrupt some disease-related chemical process. Unfortunately, drug development is an extremely long and expensive process: at the time of writing, bringing a new drug to the market typically takes more than a decade and costs more than a billion dollars. One of the reasons is the cost of testing where many drugs fail at different stages – less than 5% of candidates make it to the last stage (see e.g. Gaudelet et al. (2020)).

Since the space of chemically synthesisable molecules is very large (estimated around $10^{60}$), the search for candidate molecules with the right combination of properties such as target binding affinity, low toxicity, solubility, etc. cannot be done experimentally, and virtual or in silico screening (i.e., the use of computational techniques to identify promising molecules), is employed. Machine learning techniques play an increasingly more prominent role in this task. A prominent example of the use of Geometric Deep Learning for virtual drug screening was recently shown by Stokes et al. (2020) using a graph neural network trained to predict whether or not candidate molecules inhibit growthin the model bacterium Escherichia coli, they were able to effectively discover that Halicin, a molecule originally indicated for treating diabetes, is a highly potent antibiotic, even against bacteria strains with known antibiotic resistance. This discovery was widely covered in both scientific and popular press.

Speaking more broadly, the application of graph neural networks to molecules modeled as graphs has been a very active field, with multiple specialised architectures proposed recently that are inspired by physics and e.g. incorporate equivariance to rotations and translations (see e.g. Thomas et al. (2018); Anderson et al. (2019); Fuchs et al. (2020); Satorras et al. (2021)). Further, Bapst et al. (2020) have successfully demonstrated the utility of GNNs for predictively modelling the dynamics of glass, in a manner that outperformed the previously available physics-based models. Historically, many works in computational chemistry were precursors of modern graph neural network architectures sharing many common traits with them.

Drug Repositioning

While generating entirely novel drug candidates is a potentially viable approach, a faster and cheaper avenue for developing new therapies is drug repositioning, which seeks to evaluate already-approved drugs (either alone or in combinations) for a novel purpose. This often significantly decreases the amount of clinical evaluation that is necessary to release the drug to the market. At some level of abstraction, the action of drugs on the body biochemistry and their interactions between each other and other biomolecules can be modeled as a graph, giving rise to the concept of ‘network medicine’ coined by the prominent network scientist Albert-László Barabási and advocating the use of biological networks (such as protein-protein interactions and metabolic pathways) to develop new therapies (Barabási et al., 2011).

Geometric Deep Learning offers a modern take on this class of approaches. A prominent early example is the work of Zitnik et al. (2018), who used graph neural networks to predict side effects in a form of drug repositioning known as combinatorial therapy or polypharmacy, formulated as edge prediction in a drug-drug interaction graph. The novel coronavirus pandemic, which is largely ongoing at the time of writing this text, has sparked a particular interest in attempting to apply such approaches against COVID-19 (Gysi et al., 2020). Finally, we should note that drug repositioning is not necessarily limited to synthetic molecules: Veselkov et al. (2019) applied similar approaches to drug-like molecules contained in food (since, as we mentioned, many plant-based foods contain biological analogues of compounds used in oncological therapy). One of the authors of this text is involved in a collaboration adding a creative twist to this research, by partnering with a molecular chef that designs exciting recipes based on the ‘hyperfood’ ingredients rich in such drug-like molecules.

Protein biology

Since we have already mentioned proteins as drug targets, lets us spend a few more moments on this topic. Proteins are arguably among the most important biomolecules that have myriads of functions in our body, including protection against pathogens (antibodies), giving structure to our skin (collagen), transporting oxygen to cells (haemoglobin), catalysing chemical reactions (enzymes), and signaling (many hormones are proteins). Chemically speaking, a protein is a biopolymer, or a chain of small building blocks called aminoacids that under the influence of electrostatic forces fold into a complex 3D structure. It is this structure that endows the protein with its functions,and hence it is crucial to the understanding of how proteins work and what they do. Since proteins are common targets for drug therapies, the pharmaceutical industry has a keen interest in this field.

A typical hierarchy of problems in protein bioinformatics is going from protein sequence (a 1D string over an alphabet of of 20 different amino acids) to 3D structure (a problem known as ‘protein folding’) to function (‘protein function prediction’). Recent approaches such as DeepMind’s AlphaFold by Senior et al. (2020) used contact graphs to represent the protein structure. Gligorijevic et al. (2020) showed that applying graph neural networks on such graphs allows to achieve better function prediction than using purely sequence-based methods.

Gainza et al. (2020) developeda Geometric Deep Learning pipeline called MaSIF predicting interactions between proteins from their 3D structure. MaSIF models the protein as a molecular surface discretised as a mesh, arguing that this representation is advantageous when dealing with interactions as it allows to abstract the internal fold structure. The architecture was based on mesh convolutional neural network operating on pre-computed chemical and geometric features in small local geodesic patches. The network was trained using a few thousand co-crystal protein 3D structures from the Protein Data Bank to address multiple tasks, including interface prediction, ligand classification, and docking, and allowed to do de novo (‘from scratch’) design of proteins that could in principle act as biological immunotherapy drug against cancer – such proteins are designed to inhibit protein-protein interactions (PPI) between parts of the programmed cell death protein complex (PD-1/PD-L1) and give the immune system the ability to attack the tumor cells.

Recommender Systems and Social Networks

The first popularised large-scale applications of graph representation learning have occurred within social networks, primarily in the context of recommender systems. Recommenders are tasked with deciding which content to serve to users, potentially depending on their previous history of interactions on the service. This is typically realised through a link prediction objective: supervise the embeddings of various nodes (pieces of content) such that they are kept close together if they are deemed related (e.g. commonly viewed together). Then the proximity of two embeddings (e.g. their inner product) can be interpreted as the probability that they are linked by an edge in the content graph, and hence for any content queried by users, one approach could serve its $k$ nearest neighbours in the embedding space.

Among the pioneers of this methodology is the American image sharing and social media company Pinterest: besides presenting one of the first successful deployments of GNNs in production, their method, PinSage, successfully made graph representation learning scalable to graphs of millions of nodes and billions of edges (Ying et al., 2018). Related applications, particularly in the space of product recommendations, soon followed. Popular GNN-backed recommenders that are currently deployed in production include Alibaba’s Aligraph (Zhu et al., 2019) and Amazon’s P-Companion (Hao et al., 2020). In this way, graph deep learning is influencing millions of people on a daily level.

Within the context of content analysis on social networks, another noteworthy effort is Fabula AI, which is among the first GNN-based startups to be acquired (in 2019, by Twitter). The startup, founded by one of the authors of the text and his team, developed novel technology for detecting misinformation on social networks (Monti et al., 2019). Fabula’s solution consists of modelling the spread of a particular news item by the network of users who shared it. The users are connected if one of them re-shared the information from the other, but also if they follow each other on the social network. This graph is then fed into a graph neural network, which classifies the entire graph as either ‘true’ or ‘fake’ content – with labels based on agreement between fact-checking bodies. Besides demonstrating strong predictive power which stabilises quickly (often within a few hours of the news spreading), analysing the embeddings of individual user nodes revealed clear clustering of users who tend to share incorrect information, exemplifying the well-known ‘echo chamber’ effect.

Traffic forecasting

Transportation networks are another areawhere Geometric Deep Learning techniques are already making an actionable impact over billions of users worldwide. For example, on road networks, we can observe intersections as nodes, and road segments as edges connecting them—these edges can then be featurised by the road length, current or historical speeds along their segment, and the like.

One standard prediction problem in this space is predicting the estimated time of arrival (ETA): for a given candidate route, providing the expected travel time necessary to traverse it. Such a problem is essential in this space, not only for user-facing traffic recommendation apps, but also for enterprises (such as food delivery or ride-sharing services) that leverage these predictions within their own operations.

Graph neural networkshave shown immense promise in this space as well: they can, for example, be used to directly predict the ETA for a relevant subgraph of the road network (effectively, a graph regression task). Such an approach was successfully leveraged by DeepMind, yielding a GNN-based ETA predictor which is now deployed in production at Google Maps (Derrow-Pinion et al., 2021), serving ETA queries in several major metropolitan areas worldwide. Similar returns have been observed by the Baidu Maps team, where travel time predictions are currently served by the ConSTGAT model, which is itself based on a spatio-temporal variant of the graph attention network model (Fang et al., 2020).

Object recognition

A principal benchmark for machine learningtechniques in computer vision is the ability to classify a central object within a provided image. The ImageNet large scale visual recognition challenge (Russakovsky et al., 2015, ILSVRC) was an annual object classification challenge that propelled much of the early development in Geometric Deep Learning. ImageNet requires models to classify realistic images scraped from the Web into one of 1000 categories: such categories are at the same time diverse (covering both animate and inanimate objects), and specific (with many classes focused on distinguishing various cat and dog breeds). Hence, good performance on ImageNet often implies a solid level of feature extraction from general photographs, which formed a foundation for various transfer learning setups from pre-trained ImageNet models.

The success of convolutional neural networks on ImageNet—particularly the AlexNet model of Krizhevsky et al. (2012), which swept ILSVRC 2012 by a large margin—has in a large way spearheaded the adoption of deep learning as a whole, both in academia and in industry. Since then, CNNs have consistently ranked on top of the ILSVRC, spawning many popular architectures such as VGG-16 (Simonyan and Zisserman, 2014), Inception (Szegedy et al., 2015) and ResNets (He et al., 2016), which have successfully surpassed human-level performance on this task. The design decisions and regularisation techniques employed by these architectures (such as rectified linear activations (Nair and Hinton, 2010), dropout (Srivastava et al., 2014), skip connections (He et al., 2016) and batch normalisation (Ioffe and Szegedy, 2015)) form the backbone of many of the effective CNN models in use today.

Concurrently with object classification, significant progress had been made on object detection; that is, isolating all objects of interest within an image, and tagging them with certain classes. Such a task is relevant in a variety of downstream problems, from image captioning all the way to autonomous vehicles. It necessitates a more fine-grained approach, as the predictions need to be localised; as such, often, translation equivariant models have proven their worth in this domain. One impactful example in this space includes the R-CNN family of models (Girshick et al., 2014; Girshick, 2015; Ren et al., 2015; He et al., 2017) whereas, in the related field of semantic segmentation, the SegNet model of Badrinarayanan et al. (2017) proved influential, with its encoder-decoder architecture relying on the VGG-16 backbone.

Game playing

Convolutional neural networks also play a prominent role as translation-invariant feature extractors in reinforcement learning (RL) environments, whenever the observed state can be represented in a grid domain; e.g. this is the case when learning to play video games from pixels. In this case, the CNN is responsible for reducing the input to a flat vector representation, which is then used for deriving policy or value functions that drive the RL agent’s behaviour. While the specifics of reinforcement learning are not the focus of this section, we do note that some of the most impactful results of deep learning in the past decade have come about through CNN-backed reinforcement learning.

One particular example that is certainly worth mentioning here is DeepMind’s AlphaGo (Silver et al., 2016). It encodes the current state within a game of Go by applying a CNN to the $19\times 19$ grid representing the current positions of the placed stones. Then, through a combination of learning from previous expert moves, Monte Carlo tree search, and self-play, it had successfully reached a level of Go mastery that was sufficient to outperform Lee Sedol, one of the strongest Go players of all time, in a five-round challenge match that was widely publicised worldwide.

While this already represented a significant milestone for broader artificial intelligence—with Go having a substantially more complex state-space than, say, chess—the development of AlphaGo did not stop there. The authors gradually removed more and more Go-specific biases from the architecture, with AlphaGo Zero removing human biases, optimising purely through self-play (Silver et al., 2017), AlphaZero expands this algorithm to related two-player games, such as Chess and Shogi; lastly, MuZero (Schrittwieser et al., 2020) incorporates a model that enables learning the rules of the game on-the-fly, which allows reaching strong performance in the Atari 2600 console, as well as Go, Chess and Shogi, without any upfront knowledge of the rules. Throughout all of these developments, CNNs remained the backbone behind these models’ representation of the input.

While several high-performing RL agents were proposed for the Atari 2600 platform over the years (Mnih et al., 2015, 2016; Schulman et al., 2017), for a long time they were unable to reach human-level performance on all of the 57 games provided therein. This barrier was finally broken with Agent57 (Badia et al., 2020), which used a parametric family of policies, ranging from strongly exploratory to purely exploitative, and prioritising them in different ways during different stages of training. It, too, powers most of its computations by a CNN applied to the video game’s framebuffer.

Text and speech synthesis

Besides images (which naturally map to a two-dimensional grid), several of (geometric) deep learning’s strongest successes have happened on one-dimensional grids. Natural examples of this are text and speech, folding the Geometric Deep Learning blueprint within diverse areas such as natural language processing and digital signal processing.

Some of the most widely applied and publicised works in this space focus on synthesis: being able to generate speech or text, either unconditionally or conditioned on a particular prompt. Such a setup can support a plethora of useful tasks, such as text-to-speech (TTS), predictive text completion, and machine translation. Various neural architectures for text and speech generation have been proposed over the past decade, initially mostly based on recurrent neural networks (e.g. the aforementioned seq2seq model (Sutskever et al., 2014) or recurrent attention (Bahdanau et al., 2014)). However, in recent times, they have been gradually replaced by convolutional neural networks and Transformer-based architectures.

One particular limitation of simple 1D convolutions in this setting is their linearly growing receptive field, requiring many layers in order to cover the sequence generated so far. Dilatedconvolutions, instead, offer an exponentially growing receptive field with an equivalent number of parameters. Owing to this, they proved a very strong alternative, eventually becoming competitive with RNNs on machine translation (Kalchbrenner et al., 2016), while drastically reducing the computational complexity, owing to their parallelisability across all input positions.The most well-known application of dilated convolutions is the WaveNet model from van den Oord et al. (2016a). WaveNets demonstrated that, using dilations, it is possible to synthesise speech at the level of raw waveform (typically 16,000 samples per second or more), producing speech samples that were significantly more “human-like” than the best previous text-to-speech (TTS) systems. Subsequently, it was further demonstrated that the computations of WaveNets can be distilled in a much simpler model, the WaveRNN (Kalchbrenner et al., 2018)—and this model enabled effectively deploying this technology at an industrial scale. This allowed not only its deployment for large-scale speech generation for services such as the Google Assistant, but also allowing for efficient on-device computations; e.g. for Google Duo, which uses end-to-end encryption.

Transformers (Vaswani et al., 2017) have managed to surpass the limitations of both recurrent and convolutional architectures, showing that self-attention is sufficient for achieving state-of-the-art performance in machine translation. Subsequently, they have revolutionised natural language processing. Through the pre-trained embeddings provided by models such as BERT (Devlin et al., 2018), Transformer computations have become enabled for a large amount of downstream applications of natural language processing—for example, Google uses BERT embeddings to power its search engine.

Arguably the most widely publicised application of Transformers in the past years is text generation, spurred primarily by the Generative Pre-trained Transformer (GPT, Radford et al. (2018, 2019); Brown et al. (2020)) family of models from OpenAI. In particular, GPT-3 (Brown et al., 2020) successfully scaled language model learning to 175 billion learnable parameters, trained on next-word prediction on web-scale amounts of scraped textual corpora. This allowed it not only to become a highly-potent few-shot learner on a variety of language-based tasks, but also a text generator with capability to produce coherent and human-sounding pieces of text. This capability not only implied a large amount of downstream applications, but also induced a vast media coverage.

Healthcare

Applications in the medical domain are another promising area for Geometric Deep Learning. There are multiple ways in which these methods are being used. First, more traditional architectures such as CNNs have been applied to grid-structured data, for example, for the prediction of length of stay in Intensive Care Units (Rocheteau et al., 2020), or diagnosis of sight-threatening diseases from retinal scans (De Fauw et al., 2018). Winkels and Cohen (2019) showed that using 3D roto-translation group convolutional networks improves the accuracy of pulmonary nodule detection compared to conventional CNNs.

Second, modelling organs as geometric surfaces, mesh convolutional neural networks were shown to be able to address a diverse range of tasks, from reconstructing facial structure from genetics-related information (Mahdi et al., 2020) to brain cortex parcellation (Cucurull et al., 2018) to regressing demographic properties from cortical surface structures (Besson et al., 2020). The latter examples represent an increasing trend in neuroscience to consider the brain as a surface with complex foldsgiving rise to highly non-Euclidean structures.

At the same time, neuroscientists often try construct and analyse functional networks of the brain representing the various regions of the brain that are activated together when performing some cognitive function; these networks are often constructed using functional magnetic resonance imaging (fMRI) that shows in real time which areas of the brain consume more blood.These functional networks can reveal patient demographics (e.g., telling apart males from females, Arslan et al. (2018)), as well as used for neuropathology diagnosis, which is the third area of application of Geometric Deep Learning in medicine we would like to highlight here. In this context, Ktena et al. (2017) pioneered the use of graph neural networks for the prediction of neurological conditions such as Autism Spectrum Disorder. The geometric and functional structure of the brain appears to be intimately related, and recently Itani and Thanou (2021) pointed to the benefits of exploiting them jointly in neurological disease analysis.

Fourth, patient networks are becoming more prominent in ML-based medical diagnosis. The rationale behind these methods is that the information of patient demographic, genotypic, and phenotypic similarity could improve predicting their disease. Parisot et al. (2018) applied graph neural networks on networks of patients created from demographic features for neurological disease diagnosis, showing that the use of the graph improves prediction results. Cosmo et al. (2020) showed the benefits of latent graph learning (by which the network learns an unknown patient graph) in this setting. The latter work used data from the UK Biobank, a large-scale collection of medical data including brain imaging (Miller et al., 2016).

A wealth of data about hospital patients may be found in electronic health records (EHRs). Besides giving a comprehensive view of the patient’s progression, EHR analysis allows for relating similar patients together. This aligns with the pattern recognition method, which is commonly used in diagnostics. Therein, the clinician uses experience to recognise a pattern of clinical characteristics, and it may be the primary method used when the clinician’s experience may enable them to diagnose the condition quickly. Along these lines, several works attempt to construct a patient graph based on EHR data, either by analysing the embeddings of their doctor’s notes (Malone et al., 2018), diagnosis similarity on admission (Rocheteau et al., 2021), or even assuming a fully-connected graph (Zhu and Razavian, 2019). In all cases, promising results have been shown in favour of using graph representation learning for processing EHRs.

Particle physics and astrophysics

High energy physicists were perhaps among the first domain experts in the field of natural sciences to embrace the new shiny tool, graph neural networks. In a recent review paper, Shlomi et al. (2020)note that machine learning has historically been heavily used in particle physics experiments, either to learn complicated inverse functions allowing to infer the underlying physics process from the information measured in the detector, or to perform classification and regression tasks. For the latter, it was often necessary to force the data into an unnatural representation such as grid, in order to be able to used standard deep learning architectures such as CNN. Yet, many problems in physics involve data in the form of unordered sets with rich relations and interactions, which can be naturally represented as graphs.

One important application in high-energy physics is the reconstruction and classification of particle jets – sprays of stable particles arising from multiple successive interaction and decays of particles originating from a single initial event. In the Large Hardon Collider, the largest and best-known particle accelerator built at CERN, such jet are the result of collisions of protons at nearly the speed of light. These collisions produce massive particles, such as the long though-for Higgs boson or the top quark. The identification and classification of collision events is of crucial importance, as it might provide experimental evidence to the existence of new particles.

Multiple Geometric Deep Learning approacheshave recently been proposed for particle jet classification task, e.g. by Komiske et al. (2019) and Qu and Gouskos (2019), based on DeepSet and Dynamic Graph CNN architectures, respectively. More recently, there has also been interest in developing specialsed architectures derived from physics consideration and incorporating inductive biases consistent with Hamiltonian or Lagrangian mechanics (see e.g. Sanchez-Gonzalez et al. (2019); Cranmer et al. (2020)), equivariant to the Lorentz group (a fundamental symmetry of space and time in physics) (Bogatskiy et al., 2020), or even incorporating symbolic reasoning (Cranmer et al., 2019) and capable of learning physical laws from data. Such approaches are more interpretable (and thus considered more ‘trustworthy’ by domain experts) and also offer better generalisation.

Besides particle accelerators, particle detectors are now being used by astrophysicist for multi-messenger astronomy – a new way of coordinated observation of disparate signals, such as electromagnetic radiation, gravitational waves, and neutrinos, coming from the same source. Neutrino astronomy is of particular interest, since neutrinos interact only very rarely with matter, and thus travel enormous distances practically unaffected.Detecting neutrinos allows to observe objects inaccessible to optical telescopes, but requires enormously-sized detectors – the IceCube neutrino observatory uses a cubic kilometer of Antarctic ice shelf on the South Pole as its detector. Detecting high-energy neutrinos can possibly shed lights on some of the most mysterious objects in the Universe, such as blazars and black holes. Choma et al. (2018) used a Geometric neural network to model the irregular geometry of the IceCube neutrino detector, showing significantly better performance in detecting neutrinos coming from astrophysical sources and separating them from background events.

While neutrino astronomy offers a big promise in the study of the Cosmos, traditional optical and radio telescopes are still the ‘battle horses’ of astronomers. With these traditional instruments, Geometric Deep Learning can still offer new methodologies for data analysis. For example, Scaife and Porter (2021) used rotationally-equivariant CNNs for the classification of radio galaxies, and McEwen et al. (2021) used spherical CNNs for the analysis of cosmic microwave background radiation, a relic from the Big Bang that might shed light on the formation of the primordial Universe. As we already mentioned, such signals are naturally represented on the sphere and equivariant neural networks are an appropriate tool to study them.

Virtual and Augmented Reality

Another field of applications which served as the motivation for the development of a large class of Geometric Deep Learning methods is computer vision and graphics, in particular, dealing with 3D body models for virtual and augmented reality. Motion capture technology used to produce special effects in movies like Avatar often operates in two stages: first, the input from a 3D scanner capturing the motions of the body or the face of the actor is put into correspondence with some canonical shape, typically modelled as a discrete manifold or a mesh (this problem is often called ‘analysis’). Second, a new shape is generated to repeat the motion of the input (‘synthesis’). Initial works on Geometric Deep Learning in computer graphics and vision (Masci et al., 2015; Boscaini et al., 2016a; Monti et al., 2017) developed mesh convolutional neural networks to address the analysis problem, or more specifically, deformable shape correspondence.

First geometric autoencoder architectures for 3D shape synthesis were proposed independently by Litany et al. (2018) and Ranjan et al. (2018). In these architectures, a canonical mesh (of the body, face, or hand) was assumed to be known and the synthesis task consisted of regressing the 3D coordinates of the nodes (the embedding of the surface, using the jargon of differential geometry). Kulon et al. (2020) showed a hybrid pipeline for 3D hand pose estimation with an image CNN-based encoder and a geometric decoder. A demo of this system, developed in collaboration with a British startup company Ariel AI and presented at CVPR 2020, allowed to create realistic body avatars with fully articulated hands from video input on a mobile phone faster than real-time. Ariel AI was acquired by Snap in 2020, and at the time of writing its technology is used in Snap’s augmented reality products.

7 Historic Perspective

“Symmetry, as wide or as narrow as you may define its meaning, is one idea by which man through the ages has tried to comprehend and create order, beauty, and perfection.” This somewhat poetic definition of symmetry is given in the eponymous book of the great mathematician Hermann Weyl (2015), his Schwanengesang on the eve of retirement from the Institute for Advanced Study in Princeton. Weyl traces the special place symmetry has occupied in science and art to the ancient times, from Sumerian symmetric designs to the Pythagoreans who believed the circle to be perfect due to its rotational symmetry. Plato considered the five regular polyhedra bearing his name today so fundamental that they must be the basic building blocks shaping the material world. Yet, though Plato is credited with coining the term ςυμμετρία, which literally translates as ‘same measure’, he used it only vaguely to convey the beauty of proportion in art and harmony in music. It was the astronomer and mathematician Johannes Kepler to attempt the first rigorous analysis of the symmetric shape of water crystals. In his treatise (‘On the Six-Cornered Snowflake’),he attributed the six-fold dihedral structure of snowflakes to hexagonal packing of particles – an idea that though preceded the clear understanding of how matter is formed, still holds today as the basis of crystallography (Ball, 2011).

Symmetry in Mathematics and Physics

In modern mathematics, symmetry is almost univocally expressed in the language of group theory. The origins of this theory are usually attributed to Évariste Galois, who coined the term and used it to study solvability of polynomial equations in the 1830s. Two other names associated with group theory are those of Sophus Lie and Felix Klein, who met and worked fruitfully together for a period of time (Tobies, 2019). The former would develop the theory of continuous symmetries that today bears his name; the latter proclaimed group theory to be the organising principle of geometry in his Erlangen Program, which we mentioned in the beginning of this text. Riemannian geometry was explicitly excluded from Klein’s unified geometric picture, and it took another fifty years before it was integrated, largely thanks to the work of Élie Cartan in the 1920s.

Emmy Noether, Klein’s colleague in Göttingen, proved that every differentiable symmetry of the action of a physical system has a corresponding conservation law (Noether, 1918). In physics, it was a stunning result: beforehand, meticulous experimental observation was required to discover fundamental laws such as the conservation of energy, and even then, it was an empirical result not coming from anywhere. Noether’s Theorem — “a guiding star to 20th and 21st century physics”, in the words of the Nobel laureate Frank Wilczek — showed that the conservation of energy emerges from the translational symmetry of time, a rather intuitive idea that the results of an experiment should not depend on whether it is conducted today or tomorrow.

The symmetryassociated with charge conservation is the global gauge invariance of the electromagnetic field, first appearing in Maxwell’s formulation of electrodynamics (Maxwell, 1865); however, its importance initially remained unnoticed. The same Hermann Weyl who wrote so dithyrambically about symmetry is the one who first introduced the concept of gauge invariance in physics in the early 20th century, emphasizing its role as a principle from which electromagnetism can be derived. It took several decades until this fundamental principle — in its generalised form developed by Yang and Mills (1954) — proved successful in providing a unified framework to describe the quantum-mechanical behavior of electromagnetism and the weak and strong forces, finally culminating in the Standard Model that captures all the fundamental forces of nature but gravity. We can thus join another Nobel-winning physicist, Philip Anderson (1972), in concluding that “it is only slightly overstating the case to say that physics is the study of symmetry.”

Early Use of Symmetry in Machine Learning

In machine learning and its applications to pattern recognition and computer vision, the importance of symmetry has long been recognised. Early work on designing equivariant feature detectors for pattern recognition was done by Amari (1978),Kanatani (2012), and Lenz (1990). In the neural networks literature, the famous Group Invariance Theorem for Perceptrons by Minsky and Papert (2017) puts fundamental limitations on the capabilities of (single-layer) perceptrons to learn invariants. This was one of the primary motivations for studying multi-layer architectures (Sejnowski et al., 1986; Shawe-Taylor, 1989, 1993), which ultimately led to deep learning.

In the neural network community, Neocognitron (Fukushima and Miyake, 1982) is credited as the first implementation of shift invariance in a neural network for “pattern recognition unaffected by shift in position”. His solution came in the form of hierarchical neural network with local connectivity, drawing inspiration from the receptive fields discovered in the visual cortex by the neuroscientists David Hubel and Torsten Wiesel two decades earlier (Hubel and Wiesel, 1959). These ideas culminated in Convolutional Neural Networks in the seminal work of Yann LeCun and co-authors (LeCun et al., 1998). The first work to take a representation-theoretical view on invariant and equivariant neural networks was performed by Wood and Shawe-Taylor (1996), unfortunately rarely cited. More recent incarnations of these ideas include the works of Makadia et al. (2007); Esteves et al. (2020) and one of the authors of this text (Cohen and Welling, 2016).

Graph Neural Networks

It is difficult to pinpoint exactly when the concept of Graph Neural Networks began to emerge—partly due to the fact that most of the early work did not place graphs as a first-class citizen, partly since GNNs became practical only in the late 2010s, and partly because this field emerged from the confluence of several research areas. That being said, early forms of graph neural networks can be traced back at least to the 1990s, with examples including Alessandro Sperduti’s Labeling RAAM (Sperduti, 1994), the “backpropagation through structure” of Goller and Kuchler (1996), and adaptive processing of data structures (Sperduti and Starita, 1997; Frasconi et al., 1998). While these works were primarily concerned with operating over “structures” (often trees or directed acyclic graphs), many of the invariances preserved in their architectures are reminiscent of the GNNs more commonly in use today.

The first proper treatment of the processing of generic graph structures (and the coining of the term “graph neural network”) happened after the turn of the 21st century.Within the Artificial Intelligence lab at the Università degli Studi di Siena (Italy), papers led by Marco Gori and Franco Scarselli have proposed the first “GNN” (Gori et al., 2005; Scarselli et al., 2008). They relied on recurrent mechanisms, required the neural network parameters to specify contraction mappings, and thus computing node representations by searching for a fixed point—this in itself necessitated a special form of backpropagation (Almeida, 1990; Pineda, 1988) and did not depend on node features at all. All of the above issues were rectified by the Gated GNN (GGNN) model of Li et al. (2015). GGNNs brought many benefits of modern RNNs, such as gating mechanisms (Cho et al., 2014) and backpropagation through time, to the GNN model, and remain popular today.

Computational chemistry

It is also very important to note an independent and concurrent line of development for GNNs: one that was entirely driven by the needs of computational chemistry, where molecules are most naturally expressed as graphs of atoms (nodes) connected by chemical bonds (edges). This invited computational techniques for molecular property prediction that operate directly over such a graph structure, which had become present in machine learning in the 1990s: this includes the ChemNet model of Kireev (1995) and the work of Baskin et al. (1997). Strikingly, the “molecular graph networks” of Merkwirth and Lengauer (2005) explicitly proposed many of the elements commonly found in contemporary GNNs—such as edge type-conditioned weights or global pooling—as early as 2005. The chemical motivation continued to drive GNN development into the 2010s, with two significant GNN advancements centered around improving molecular fingerprinting (Duvenaud et al., 2015) and predicting quantum-chemical properties (Gilmer et al., 2017) from small molecules. At the time of writing this text, molecular property prediction is one of the most successful applications of GNNs, with impactful results in virtual screening of new antibiotic drugs (Stokes et al., 2020).

Node embeddings

Some of the earliest success stories of deep learning on graphs involve learning representations of nodes in an unsupervised fashion, based on the graph structure. Given their structural inspiration, this direction also provides one of the most direct links between graph representation learning and network science communities. The key early approaches in this space relied on random walk-based embeddings: learning node representations in a way that brings them closer together if the nodes co-occur in a short random walk. Representative methods in this space include DeepWalk (Perozzi et al., 2014), node2vec (Grover and Leskovec, 2016) and LINE (Tang et al., 2015), which are all purely self-supervised. Planetoid (Yang et al., 2016) was the first in this space to incorporate supervision label information, when it is available.

Unifying random walk objectives with GNN encoderswas attempted on several occasions, with representative approaches including Variational Graph Autoencoder (VGAE, Kipf and Welling (2016b)), embedding propagation (García-Durán and Niepert, 2017), and unsupervised variants of GraphSAGE (Hamilton et al., 2017). However, this was met with mixed results, and it was shortly discovered that pushing neighbouring node representations together is already a key part of GNNs’ inductive bias. Indeed, it was shown that an untrained GNN was already showing performance that is competitive with DeepWalk, in settings where node features are available (Veličković et al., 2019; Wu et al., 2019). This launched a direction that moves away from combining random walk objectives with GNNs and shifting towards contrastive approaches inspired by mutual information maximisation and aligning to successful methods in the image domain. Prominent examples of this direction include Deep Graph Informax (DGI, Veličković et al. (2019)), GRACE (Zhu et al., 2020), BERT-like objectives (Hu et al., 2020) and BGRL (Thakoor et al., 2021).

Probabilistic graphical models

Graph neural networks have also, concurrently, resurged through embedding the computations of probabilistic graphical models (PGMs, Wainwright and Jordan (2008)). PGMs are a powerful tool for processing graphical data, and their utility arises from their probabilistic perspective on the graph’s edges: namely, the nodes are treated as random variables, while the graph structure encodes conditional independence assumptions, allowing for significantly simplifying the calculation and sampling from the joint distribution. Indeed, many algorithms for (exactly or approximately) supporting learning and inference on PGMs rely on forms of passing messages over their edges (Pearl, 2014), with examples including variational mean-field inference and loopy belief propagation (Yedidia et al., 2001; Murphy et al., 2013).

This connection between PGMs and message passing was subsequently developed into GNN architectures, with early theoretical links established by the authors of structure2vec (Dai et al., 2016). Namely, by posing a graph representation learning setting as a Markov random field (of nodes corresponding to input features and latent representations), the authors directly align the computation of both mean-field inference and loopy belief propagation to a model not unlike the GNNs commonly in use today.

The key “trick” which allowed for relating the latent representations of a GNN to probability distributions maintained by a PGM was the usage of Hilbert-space embeddings of distributions (Smola et al., 2007). Given $\phi$, an appropriately chosen embedding function for features $\mathbf{x}$, it is possible to embed their probability distribution $p(\mathbf{x})$ as the expected embedding $\mathbb{E}_{\mathbf{x}\sim p(\mathbf{x})}\phi(\mathbf{x})$. Such a correspondence allows us to perform GNN-like computations, knowing that the representations computed by the GNN will always correspond to an embedding of some probability distribution over the node features.

The structure2vec model itself is, ultimately, a GNN architecture which easily sits within our framework, but its setup has inspired a series of GNN architectures which more directly incorporate computations found in PGMs. Emerging examples have successfully combined GNNs with conditional random fields (Gao et al., 2019; Spalević et al., 2020), relational Markov networks (Qu et al., 2019) and Markov logic networks (Zhang et al., 2020).

The Weisfeiler-Lehman formalism

The resurgence of graph neural networks was followed closely by a drive to understand their fundamental limitations, especially in terms of expressive power. While it was becoming evident that GNNs are a strong modelling tool of graph-structured data, it was also clear that they wouldn’t be able to solve any task specified on a graph perfectly.A canonical illustrative example of this is deciding graph isomorphism: is our GNN able to attach different representations to two given non-isomorphic graphs? This is a useful framework for two reasons. If the GNN is unable to do this, then it will be hopeless on any task requiring the discrimination of these two graphs. Further, it is currently not known if deciding graph isomorphism is in P, the complexity class in which all GNN computations typically reside.

The key framework which binds GNNs to graph isomorphism is the Weisfeiler-Lehman (WL) graph isomorphism test (Weisfeiler and Leman, 1968). This test generates a graph representation by iteratively passing node features along the edges of the graph, then randomly hashing their sums across neighbourhoods. Connections to randomly-initialised convolutional GNNs are apparent, and have been observed early on: for example, within the GCN model of Kipf and Welling (2016a). Aside from this connection, the WL iteration was previously introduced in the domain of graph kernels by Shervashidze et al. (2011), and it still presents a strong baseline for unsupervised learning of whole-graph representations.

Whilethe WL test is conceptually simple, and there are many simple examples of non-isomorphic graphs it cannot distinguish, its expressive power is ultimately strongly tied to GNNs. Analyses by Morris et al. (2019) and Xu et al. (2018) have both reached a striking conclusion: any GNN conforming to one of the three flavours we outlined in Section 5.3 cannot be more powerful than the WL test!

In order to exactly reach this level of representational power, certain constraints must exist on the GNN update rule. Xu et al. (2018) have shown that, in the discrete-feature domain, the aggregation function the GNN uses must be injective, with summation being a key representative. Based on the outcome of their analysis, Xu et al. (2018) propose the Graph Isomorphism Network (GIN), which is a simple but powerful example of a maximally-expressive GNN under this framework. It is also expressible under the convolutional GNN flavour we propose.

Lastly, it is worth noting that these findings do not generalise to continuous node feature spaces. In fact, using the Borsuk-Ulam theorem (Borsuk, 1933), Corso et al. (2020) have demonstrated that, assuming real-valued node features, obtaining injective aggregation functions requires multiple aggregators (specifically, equal to the degree of the receiver node). Their findings have driven the Principal Neighbourhood Aggregation (PNA) architecture, which proposes a multiple-aggregator GNN that is empirically powerful and stable.

Higher-order methods

The findings of the previous paragraphs do not contradict the practical utility of GNNs. Indeed, in many real-world applications the input features are sufficiently rich to support useful discriminative computations over the graph structure, despite of the above limitations.

However, one key corollary is that GNNs are relatively quite weak at detecting some rudimentary structures within a graph. Guided by the specific limitations or failure cases of the WL test, several works have provided stronger variants of GNNs that are provably more powerful than the WL test, and hence likely to be useful on tasks that require such structural detection.

Perhaps the most direct place to hunt for more expressive GNNs is the WL test itself. Indeed, the strength of the original WL test can be enhanced by considering a hierarchy of WL tests, such that $k$-WL tests attach representations to $k$-tuples of nodes (Morris et al., 2017). The $k$-WL test has been directly translated into a higher-order $k$-GNN architecture by Morris et al. (2019),which is provably more powerful than the GNN flavours we considered before. However, its requirement to maintain tuple representations implies that, in practice, it is hard to scale beyond $k=3$.

Concurrently, Maron et al. (2018, 2019) have studied the characterisation of invariant and equivariant graph networks over $k$-tuples of nodes. Besides demonstrating the surprising result of any invariant or equivariant graph network being expressible as a linear combination of a finite number of generators—the amount of which only depends on $k$—the authors showed that the expressive power of such layers is equivalent to the $k$-WL test, and proposed an empirically scalable variant which is provably 3-WL powerful.

Besides generalising the domain over which representations are computed, significant effort had also went into analysing specific failure cases of 1-WL and augmenting GNN inputs to help them distinguish such cases. One common example is attaching identifying features to the nodes, which can help detecting structure. Proposals to do this include one-hot representations (Murphy et al., 2019), as well as purely random features (Sato et al., 2020).

More broadly, there have been many efforts to incorporate structural information within the message passing process, either by modulating the message function or the graph that the computations are carried over. Several interesting lines of work here involve sampling anchor node sets (You et al., 2019), aggregating based on Laplacian eigenvectors (Stachenfeld et al., 2020; Beaini et al., 2020; Dwivedi and Bresson, 2020), or performing topological data analysis, either for positional embeddings (Bouritsas et al., 2020) or driving message passing (Bodnar et al., 2021).

Signal processing and Harmonic analysis

Since the early successes of Convolutional Neural Networks, researchers have resorted to tools from harmonic analysis, image processing, and computational neuroscience trying to provide a theoretical framework that explains their efficiency. $M$-theory is a framework inspired by the visual cortex, pioneered by Tomaso Poggio and collaborators (Riesenhuber and Poggio, 1999; Serre et al., 2007), based on the notion of templates that can be manipulated under certain symmetry groups. Another notable model arising from computational neuroscience were steerable pyramids, a form of multiscale wavelet decompositions with favorable properties against certain input transformations, developed by Simoncelli and Freeman (1995). They were a central element in early generative models for textures (Portilla and Simoncelli, 2000), which were subsequently improved by replacing steerable wavelet features with deep CNN features Gatys et al. (2015). Finally, Scattering transforms, introduced by Stéphane Mallat (2012) and developed by Bruna and Mallat (2013), provided a framework to understand CNNs by replacing trainable filters with multiscale wavelet decompositions, also showcasing the deformation stability and the role of depth in the architecture.

Signal Processing on Graph and Meshes

Another important class of graph neural networks, often referred to as spectral, has emerged from the work of one of the authors of this text (Bruna et al., 2013), using the notion of the Graph Fourier transform. The roots of this construction are in the signal processing and computational harmonic analysis communities, where dealing with non-Euclidean signals has become prominent in the late 2000s and early 2010s. Influential papers from the groups of Pierre Vandergheynst (Shuman et al., 2013) and José Moura (Sandryhaila and Moura, 2013) popularised the notion of “Graph Signal Processing” (GSP) and the generalisation of Fourier transforms based on the eigenvectors of graph adjacency and Laplacian matrices. The graph convolutional neural networks relying on spectral filters by Defferrard et al. (2016) and Kipf and Welling (2016a) are among the most cited in the field and can likely be credited) as ones reigniting the interest in machine learning on graphs in recent years.

It is worth noting that, in the field of computer graphics and geometry processing, non-Euclidean harmonic analysis predates Graph Signal Processing by at least a decade. We can trace spectral filters on manifolds and meshes to the works of Taubin et al. (1996). These methods became mainstream in the 2000s following the influential papers of Karni and Gotsman (2000) on spectral geometry compression and of Lévy (2006) on using the Laplacian eigenvectors as a non-Euclidean Fourier basis. Spectral methods have been used for a range of applications,most prominent of which is the construction of shape descriptors (Sun et al., 2009) and functional maps (Ovsjanikov et al., 2012); these methods are still broadly used in computer graphics at the time of writing.

Computer Graphics and Geometry Processing

Models for shape analysis based on intrinsic metric invariants were introduced by various authors in the field of computer graphics and geometry processing (Elad and Kimmel, 2003; Mémoli and Sapiro, 2005; Bronstein et al., 2006), and are discussed in depth by one of the authors in his earlier book (Bronstein et al., 2008). The notions of intrinsic symmetries were also explored in the same field Raviv et al. (2007); Ovsjanikov et al. (2008). The first architecture for deep learning on meshes, Geodesic CNNs, was developed in the team of one of the authors of the text (Masci et al., 2015). This model used local filters with shared weights, applied to geodesic radial patches. It was a particular setting of gauge-equivariant CNNs developed later by another author of the text (Cohen et al., 2019). A generalisation of Geodesic CNNs with learnable aggregation operations, MoNet, proposed by Federico Monti et al. (2017) from the same team, used an attention-like mechanism over the local structural features of the mesh, that was demonstrated to work on general graphs as well. The graph attention network (GAT), which technically speaking can be considered a particular instance of MoNet, was introduced by another author of this text (Veličković et al., 2018). GATs generalise MoNet’s attention mechanism to also incorporate node feature information, breaking away from the purely structure-derived relevance of prior work. It is one of the most popular GNN architectures currently in use.

In the context of computer graphics, it is also worthwhile to mention that the idea of learning on sets (Zaheer et al., 2017) was concurrently developed in the group of Leo Guibas at Stanford under the name PointNet (Qi et al., 2017) for the analysis of 3D point clouds. This architecture has lead to multiple follow-up works, including one by an author of this text called Dynamic Graph CNN (DGCNN, Wang et al. (2019b)). DGCNN used a nearest-neighbour graph to capture the local structure of the point cloud to allow exchange of information across the nodes; the key characteristic of this architecture was that the graph was constructed on-the-fly and updated between the layers of the neural network in relation to the downstream task. This latter property made DGCNN one of the first incarnations of ‘latent graph learning’, which in its turn has had significant follow up. Extensions to DGCNN’s $k$-nearest neighbour graph proposal include more explicit control over these graphs’ edges, either through bilevel optimisation (Franceschi et al., 2019), reinforcement learning (Kazi et al., 2020) or direct supervision (Veličković et al., 2020). Independently, a variational direction (which probabilistically samples edges from a computed posterior distribution) has emerged through the NRI model (Kipf et al., 2018). While it still relies on quadratic computation in the number of nodes, it allows for explicitly encoding uncertainty about the chosen edges.

Another very popular direction in learning on graphs without a provided graph relies on performing GNN-style computations over a complete graph, letting the network infer its own way to exploit the connectivity. The need for this arisen particularly in natural language processing, where various words in a sentence interact in highly nontrivial and non-sequential ways. Operating over a complete graph of words brought about the first incarnation of the Transformer model (Vaswani et al., 2017), which de-throned both recurrent and convolutional models as state-of-the-art in neural machine translation, and kicked off an avalanche of related work, transcending the boundaries between NLP and other fields. Fully-connected GNN computation has also concurrently emerged on simulation (Battaglia et al., 2016), reasoning (Santoro et al., 2017), and multi-agent (Hoshen, 2017) applications, and still represents a popular choice when the number of nodes is reasonably small.

Algorithmic reasoning

For most of the discussion we posed in this section, we have given examples of spatially induced geometries, which in turn shape the underlying domain, and its invariances and symmetries. However, plentiful examples of invariances and symmetries also arise in a computational setting. One critical difference to many common settings of Geometric Deep Learning is that links no longer need to encode for any kind of similarity, proximity, or types of relations—they merely specify the “recipe” for the dataflow between data points they connect.

Instead, the computations of the neural network mimic the reasoning process of an algorithm (Cormen et al., 2009), with additional invariances induced by the algorithm’s control flow and intermediate results. In the space of algorithms the assumed input invariants are often referred to as preconditions, while the invariants preserved by the algorithm are known as postconditions.

Eponymously, the research direction of algorithmic reasoning (Cappart et al., 2021, Section 3.3.) seeks to produce neural network architectures that appropriately preserve algorithmic invariants. The area has investigated the construction of general-purpose neural computers, e.g., the neural Turing machine (Graves et al., 2014) and the differentiable neural computer (Graves et al., 2016). While such architectures have all the hallmarks of general computation, they introduced several components at once, making them often challenging to optimise, and in practice, they are almost always outperformed by simple relational reasoners, such as the ones proposed by Santoro et al. (2017, 2018).

As modelling complex postconditions is challenging, plentiful work on inductive biases for learning to execute (Zaremba and Sutskever, 2014) has focused on primitive algorithms (e.g. simple arithmetic). Prominent examples in this space include the neural GPU (Kaiser and Sutskever, 2015), neural RAM (Kurach et al., 2015), neural programmer-interpreters (Reed and De Freitas, 2015), neural arithmetic-logic units (Trask et al., 2018; Madsen and Johansen, 2020) and neural execution engines (Yan et al., 2020).

Emulating combinatorial algorithms of superlinear complexity was made possible with the rapid development of GNN architectures. The algorithmic alignment framework pioneered by Xu et al. (2019) demonstrated, theoretically, that GNNs align with dynamic programming (Bellman, 1966), which is a language in which most algorithms can be expressed. It was concurrently empirically shown, by one of the authors of this text, that it is possible to design and train GNNs that align with algorithmic invariants in practice (Veličković et al., 2019). Onwards, alignment was achieved with iterative algorithms (Tang et al., 2020), linearithmic algorithms (Freivalds et al., 2019), data structures (Veličković et al., 2020) and persistent memory (Strathmann et al., 2021). Such models have also seen practical use in implicit planners (Deac et al., 2020), breaking into the space of reinforcement learning algorithms.

Concurrently, significant progress has been made on using GNNs for physics simulations (Sanchez-Gonzalez et al., 2020; Pfaff et al., 2020). This direction yielded much of the same recommendations for the design of generalising GNNs. Such a correspondence is to be expected: given that algorithms can be phrased as discrete-time simulations, and simulations are typically implemented as step-wise algorithms, both directions will need to preserve similar kinds of invariants.

Tightly bound with the study of algorithmic reasoning are measures of extrapolation. This is a notorious pain-point for neural networks, given that most of their success stories are obtained when generalising in-distribution; i.e. when the patterns found in the training data properly anticipate the ones found in the test data. However, algorithmic invariants must be preserved irrespective of, e.g., the size or generative distribution of the input, meaning that the training set will likely not cover any possible scenario encountered in practice. Xu et al. (2020b) have proposed a geometric argument for what is required of an extrapolating GNN backed by rectifier activations: its components and featurisation would need to be designed so as to make its constituent modules (e.g. message function) learn only linear target functions. Bevilacqua et al. (2021) propose observing extrapolation under the lens of causal reasoning, yielding environment-invariant representations of graphs.

Geometric Deep Learning

Our final historical remarks regard the very name of this text. The term ‘Geometric Deep Learning’ was first introduced by one of the authors of this text in his ERC grant in 2015 and popularised in the eponymous IEEE Signal Processing Magazine paper (Bronstein et al., 2017). This paper proclaimed, albeit “with some caution”, the signs of “a new field being born.” Given the recent popularity of graph neural networks, the increasing use of ideas of invariance and equivariance in a broad range of machine learning applications, and the very fact of us writing this text, it is probably right to consider this prophecy at least partially fulfilled. The name “4G: Grids, Graphs, Groups, and Gauges” was coined by Max Welling for the ELLIS Program on Geometric Deep Learning, co-directed by two authors of the text. Admittedly, the last ‘G’ is somewhat of a stretch, since the underlying structures are manifolds and bundles rather than gauges. For this text, we added another ‘G’, Geodesics, in reference to metric invariants and intrinsic symmetries of manifolds.